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Abstract. We describe the Einstein Toolkit, a community-driven, freely accessible 
computational infrastructure intended for use in numerical relativity, relativistic 
astrophysics, and other applications. The Toolkit, developed by a collaboration involving 
researchers from multiple institutions around the world, combines a core set of components 
needed to simulate astrophysical objects such as black holes, compact objects, and collapsing 
stars, as well as a full suite of analysis tools. The Einstein Toolkit is currently based 
on the Cactus Framework for high-performance computing and the Carpet adaptive mesh 
refinement driver. It implements spacetime evolution via the BSSN evolution system and 
general-relativistic hydrodynamics in a finite-volume discretization. The toolkit is under 
continuous development and contains many new code components that have been publicly 
released for the first time and are described in this article. We discuss the motivation behind 
the release of the toolkit, the philosophy underlying its development, and the goals of the 
project. A summary of the implemented numerical techniques is included, as are results of 
numerical test covering a variety of sample astrophysical problems. 
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1. Introduction 

Scientific progress in the field of numerical relativity has always been closely tied to the 
availability and ease-of-use of enabling software and computational infrastructure. This 
document describes the Einstein Toolkit, which provides such an infrastructure, developed 
openly and made available freely with grant support from the National Science Foundation. 

Now is a particularly exciting time for numerical relativity and relativistic astrophysics, 
with major advances having been achieved in the study of astrophysical systems containing 
black holes (BHs) and neutron stars (NSs). The first fully general relativistic (GR) 
simulations of merging NS-NS binaries were reported in 1999, with further advances over 
the next few years P-[5]. However, systems containing BHs proved much more difficult to 
evolve numerically until 2005. That year, computational breakthroughs were made following 
the development of a generalized harmonic formulation |£j and then a "moving puncture" 
approach [71 18] in the BSSN (Baumgarte-Shapiro-Shibata-Nakamura) formalism [9l |10] that 
lead to the first stable long-term evolutions of moving single and multiple BH systems. These 
results quickly transformed the field which was now able to effectively evolve the Einstein 
field equations for coalescing BH-BH binaries and other systems containing moving BHs, 
including merging BH-NS binaries. 

These breakthroughs had direct relevance to astrophysics, and enabled exciting 
new results on recoil velocities from BH-BH mergers (e.g, piTHTB] and references 
therein), post-Newtonian (PN) and numerical waveform comparisons and waveform 
template generation (e.g., [TTH^ and references therein), comparisons between numerical 
waveforms [261 EZ], determination of the spin of the remnant BH formed in BH-BH 
mergers (e.g, [281433] and references therein), and studies of eccentric BH-BH binaries [3^39] . 

Meanwhile, general relativistic magneto-hydrodynamics (GRMHD) on fixed background 
spacetimes has been successful in multi-dimensional settings since the mid-1990s, focusing 
on BH accretion processes and relativistic jet production and evolution (see ^D] for a review 
of the numerical formalism and [H] for a review of work on disk and jet models). GRMHD 
coupled with curvature evolution, on the other hand, which is crucial for modeling large- 
scale bulk dynamics in compact binary star coalescence or single-star collapse scenarios, 
has started to produce astrophysically interesting results only in the past ~ 3 — 5 years, 
enabled primarily by the availability of long-term stable curvature evolution systems as well 
as improved GRMHD algorithms (see [ID] for a review). In addition to these developments, 
substantial progress has been made in importing more physically motivated equations of 
state (EOS), including tabulated versions (e.g., [l2] - [ii] ) and temperature-dependent models 
(e.g., [I5lfl7] ). Some codes have also begun to incorporate microphysical effects of neutrino 
emission and deleptonization H9] . 

Many of the successful techniques used to evolve BH-BH binaries have proven to be 
equally applicable to merging NS-NS [IHHSfflSB] and BH-NS [SHIESHHI] binaries (for reviews, 
see also [85l 186]). allowing for further investigations into the former and the first full GR 
simulations of the latter. All recent results use either the general harmonic formalism or the 
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BSSN formalism in the "moving puncture" gauge. Nearly all include some form of adaptive 
mesh refinement, since unigrid models cannot produce accurate long-term evolutions without 
requiring exorbitant computational resources, though some BH-NS simulations have been 
performed with a pseudospectral code [TOl [TU [TH [75]. Many groups' codes now include 
GRMHD (used widely for NS-NS mergers, and for BH-NS mergers in [02], and some include 
microphysical effects as well (e.g., [i8l 1661 ITT] ). 

In addition to studying binary mergers, numerical relativity is a necessary element for 
understanding stellar collapse and dynamical instabilities in NSs. GRHD has been used to 
study, among many other applications, massive stars collapsing to protoneutron stars [871 ^ 
EH], the collapse of rotating, hypermassive NSs to BHs in 2D and 3D (see, e.g., PUH55] ). and 
non-axisymmetric instabilities in rapidly rotating polytropic NS models [9T | l96 | 197]. 

Simultaneously with the advances in both our physical understanding of relativistic 
dynamics and the numerical techniques required to study them, a set of general computational 
tools and libraries has been developed with the aim of providing a computational core that 
can enable new science, broaden the community, facilitate collaborative and interdisciplinary 
research, promote software reuse and take advantage of emerging petascale computers 
and advanced cyberinfrastructure: the Cactus computational toolkit [98]. Although the 
development of Cactus was driven directly from the numerical relativity community, it 
was developed in collaboration with computer scientists and other computational fields to 
facilitate the incorporation of innovations in computational theory and technology. 

This success prompted usage of the Cactus computational toolkit in other areas, such 
as ocean forecast models [99j and chemical reaction simulations |lUUj . At the same time, 
the growing number of results in numerical relativity increased the need for commonly 
available utilities such as comparison and analysis tools, typically those specifically designed 
for astrophysical problems. Including them within the Cactus computational toolkit was 
not felt to fit within its rapidly expanding scope. This triggered the creation of the Einstein 
Toolkit |lUlj . Large parts of the Einstein toolkit presently do make use of the Cactus toolkit, 
but this is not an requirement, and other contributions are welcome, encouraged and have 
been accepted in the past. 

2. Requirements 

2.1. Scientific 

While the aforementioned studies collectively represent breakthrough simulations that have 
significantly advanced the modeling of relativistic astrophysical systems, all simulations are 
presently missing one or more critical physical ingredients and are lacking the numerical 
precision to accurately and realistically model the large-scale and small-scale dynamics of 
their target systems simultaneously. 

One of the aims of the Einstein Toolkit is to provide or extend some of these missing 
ingredients in the course of its development. Over the past three years, routines have been 
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added to the code to allow for a wider range of initial data choices, to allow for multithreading 
in hydrodynamic evolutions, and to refine the Carpet adaptive mesh refinement driver. 
Looking forward, three possible additions to future releases are the inclusion of magnetic 
fields into the dynamics via an ideal MHD treatment, more physical nuclear matter equations 
of state (EOSs) including the ability to model finite-temperature effects, and higher-order 
numerical techniques. All of these are under active development, with MHD and finite- 
temperature evolution code already available, though not completely documented, within 
the pubhc toolkit releases, and will be made available once they are thoroughly tested and 
validated against known results. 

2.2. Academic and Social 

A primary concern for research groups is securing reliable funding to support graduate 
students and postdoctoral researchers. This is easier to achieve if it can be shown that 
scientific goals can be attacked directly with fewer potential infrastructure problems, one of 
the goals of the Einstein Toolkit. 

While the Einstein Toolkit does have a large group of users, many of them do not 
directly collaborate on science problems, and some compete. However, many groups agree 
that sharing the development of the underlying infrastructure is mutually beneficial for every 
group and the wider community as well. This is achieved by lifting off the research groups' 
shoulders much of the otherwise necessary burden of creating such an infrastructure, while at 
the same time increasing the amount of code review and thus, code quality. In addition, the 
Einstein Toolkit provides computer scientists an ideal platform to perform state-of-the-art 
research, which directly benefits research in other areas of science and provides an immediate 
application of their research. 

3. Design and Strategy 

The mechanisms for the development and support of the Einstein Toolkit are designed to 
be open, transparent and community-driven. The complete source code, documentation 
and tools included in the Einstein Toolkit are distributed under open-source licenses. The 
Einstein Toolkit maintains a version control system (svn . einsteintoolkit . org) with open 
access that contains software supported by the Einstein Toolkit, the toolkit web pages, and 
documentation. An open wiki for documentation (docs.einsteintoollkit.org) has been 
established where the community can contribute either anonymously or through personal 
authentication. Almost all discussions about the toolkit take place on an open mail list 
(users@einsteintoolkit.org). The regular weekly meetings for the Einstein Toolkit are 
open and the community is invited to participate. Meeting minutes are recorded and publicly 
available as well. The Einstein Toolkit blog requires users to first request a login, but then 
allows for posting at will. Any user can post comments to entries already on the blog. The 
community makes heavy use of an issue tracking system (trac . einsteintoolkit . org), with 



The Einstein Toolkit 



5 



submissions also open to the public. 

Despite this open design, some actions naturally have to be restricted to a smaller group 
of maintainers. This is true for administrative tasks like the setup and maintenance of the 
services themselves, or to avoid large amounts of spam. One of the most important tasks 
of an Einstein Toolkit maintainer is to review and apply patches sent by users in order to 
ensure a high software quality level. Every substantial change or addition to the toolkit must 
be reviewed by another Einstein Toolkit maintainer, and is generally open for discussion on 
the users mailing list. This convention, though not being technically enforced, works well in 
practice and promotes active development. 

4. Core Technologies 

The Einstein Toolkit modules center around a set of core modules that provide basic 
functionality to create, deploy and manage a numerical simulation starting with code 
generation all to way to archiving of simulation results: (i) the Cactus framework "flesh" 
provides the underlying infrastructure to build complex simulation codes out of independently 
developed modules and facilities communication between these modules, (ii) the adaptive 
mesh refinement driver. Carpet, is build on top of Cactus and provides problem independent 
adaptive mesh refinement support for simulations that need to resolve physics on length scales 
differing by many orders of magnitude, while relieving the scientist of the need to worry about 
internal details of the mesh refinement driver, (iii) Kranc, which generates code in a computer 
language from a high-level description in Mathematica and (iv) the Simulation Factory, which 
provides a uniform, high-level interface to common operations, such as submission and restart 
of jobs, for a large number of compute clusters. 

4-1- Cactus Framework 

The Cactus Framework [HSl HU^l I1U3] is an open source, modular, portable programming 
environment for collaborative HPC computing primarily developed at Louisiana State 
University, which originated at the Albert Einstein Institute and also has roots at the National 
Center for Supercomputing Applications (see, e.g., jl04H106j for historical reviews). The 
Cactus computational toolkit consists of general modules which provide parallel drivers, 
coordinates, boundary conditions, interpolators, reduction operators, and efficient I/O in 
different data formats. Generic interfaces make it possible to use external packages and 
improved modules which are made immediately available to users. 

The structure of the Cactus framework is completely modular, with only a very small 
core (the "flesh") which provides the interfaces between modules both at compile- and 
run-time. The Cactus modules (called "thorns") may (and typically do) specify inter- 
module dependencies, e.g., to share or extend conflguration information, common variables, 
or runtime parameters. Modules compiled into an executable can remain dormant at run- 
time. This usage of modules and a common interface between them enables researchers to 
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1) easily use modules written by others without the need to understand all details of their 
implementation and 2) write their own modules without the need to change the source code 
of other parts of a simulation in the (supported) programming language of their choice. The 
number of active modules within a typical Cactus simulation ranges from tens to hundreds 
and often has an extensive set of inter-module dependencies. 

The Cactus Framework was developed originally by the numerical relativity community, 
and although it is now a general component framework that supports different application 
domains, its core user group continues to be comprised of numerical relativists. It is not 
surprising therefore, that one of the science modules provided in the Einstein Toolkit is a set 
of state of the art modules to simulate binary black hole mergers. All modules to simulate 
and analyze the data are provided out of the box. This set of modules also provides a 
way of testing the Einstein Toolkit modules in a production type simulation rather than 
synthetic test cases. Some of these modules have been developed specifically for the Einstein 
Toolkit while others are modules used in previous publications and have been contributed 
to the toolkit. In these cases the Einstein Toolkit provides documentation and best practice 
guidelines for the contributed modules. 

4-2. Adaptive Mesh Refinement 

In Cactus, infrastructure capabilities such as memory management, parallelization, time 
evolution, mesh refinement, and I/O are delegated to a set of special driver components. This 
helps separate physics code from infrastructure code; in fact, a typical physics component 
(implementing, e.g., the Einstein or relativistic MHD equations) does not contain any code 
or subroutine calls having to do with parallelization, time evolution, or mesh refinement. 
The information provided in the interface declarations of the individual components allows 
a highly efficient execution of the combined program. 

The Einstein Toolkit offers two drivers, PUGH and Carpet. PUGH provides domains 
consisting of a uniform grid with Cartesian topology, and is highly scalable (up to more 
than 130,000 |107j .) Carpet |108H110] provides multi-block methods and adaptive mesh 
refinement (AMR). Multi-block methods cover the domain with a set of (possibly distorted) 
blocks that exchange boundary information via techniques such as interpolation or penalty 
methods. I The AMR capabilities employ the standard Berger-Oliger algorithm with 
subcycling in time. 

AMR implies that resolution in the simulation domain is dynamically adapted to the 
current state of the simulation, i.e., regions that require a higher resolution are covered with 
blocks with a finer grid (typically by a factor of two); these are called refined levels. Finer 
grids can be also recursively refined. At regular intervals, the resolution requirements in the 
simulation are re-evaluated, and the grid hierarchy is updated; this step is called regridding. 

Since a finer grid spacing also requires smaller time steps for hyperbolic problems, the 

I Although muhi-block methods are supported by Carpet, the Einstein Toolkit itself does not currently 
contain any multi-block coordinate systems. 
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Figure 1. Berger-Oliger time stepping details, showing a coarse and a fine grid, with time 
advancing upwards. Left: Time stepping algorithm. First the coarse grid takes a large 
time step, then the refined grid takes two smaller steps. The fine grid solution is then 
injected into the coarse grid where the grids overlap. Right: Fine grid boundary conditions. 
The boundary points of the refined grids are filled via interpolation. This may require 
interpolation in space and in time. 



' — \ — • — • — • — • 
/M ■ ■ ■ ■ 



finer grids perform multiple time steps for eacli coarse grid time step, leading to a recursive 
time evolution pattern that is typical for Berger-Oliger AMR. If a simulation uses 11 levels, 
then the resolutions (both in space and time) of the the coarsest and finest levels differ by a 
factor of 2^^^^ = 1024. This non-uniform time stepping leads to a certain complexity that is 
also handled by the Carpet driver; for example, applying boundary conditions to a fine level 
requires interpolation in space and time from a coarser level. Outputting the solution at a 
time in between coarse grid time steps also requires interpolation. These parallel interpolation 
operations are implemented efficiently in Carpet and are applied automatically as specified in 
the execution schedule, i.e. without requiring function calls in user code. Figure [T] describes 
some details of the Berger-Oliger time stepping algorithm; more details are described in |108] . 

Carpet is the main driver used today for Cactus-based astrophysical simulations. Carpet 
offers hybrid MPI/OpenMP parallelization and is used in production on up to several 
thousand cores |112t 1113] . Figure [2] shows a weak scaling test of Carpet, where McLachlan 
(see section 5.3 below) solves the Einstein equations on a grid structure with nine levels of 



mesh refinement. This demonstrates excellent scalability up to more than ten thousand cores. 
In production simulations, smaller and more complex grid structures, serial tasks related to 
online data analysis and other necessary tasks reduce scalability by up to a factor of ten. 

We estimate that in 2010, about 7,000 core years of computing time (45 million core 
hours) will have been used via Carpet by more than a dozen research groups world-wide. To 
date, more than 90 peer-reviewed publications and more than 15 student theses have been 
based on Carpet |110] . 



4-3. Simulation Factory 

Today's supercomputers differ significantly in their hardware configuration, available 
software, directory structure, queuing system, queuing policy, and many other user- 
visible properties. In addition, the system architectures and user interfaces offered by 
supercomputers are very different from those offered by laptops or workstations. This 
makes performing large, three-dimensional time-dependent simulations a complex, time- 
consuming and difficult task. However, most of these differences are only superficial, and the 
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Figure 2. Results from weak scaling tests evolving the Einstein equations on a mesh 
refinement grid structure with nine levels. This shows the time required per grid point, where 
smaller numbers are better (the ideal scaling is a horizontal line) . This demonstrates excellent 
scalability to up to more than 10,000 cores. Including hydrodynamics approximately doubles 
calculation times without negatively influencing scalability. 



basic capabilities of supercomputers are very similar; most of the complexity of managing 
simulations lies in menial tasks that require no physical or numerical insight. 

The Simulation Factory |114t 1115] offers a set of abstractions for the tasks necessary 
to set up and successfully complete numerical simulations based on the Cactus framework. 
These abstractions hide tedious low-level management operations, capture "best practices" 
of experienced users, and create a log trail ensuring repeatable and well-documented scientific 
results. Using these abstractions, most operations are simplified and many types of 
potentially disastrous user errors are avoided, allowing different supercomputers to be used 
in a uniform manner. 

Using the Simulation Factory, we offer a tutorial for the Einstein Toolkit |lUlj that 
teaches new users how to download, configure, build, and run full simulations of the 
coupled Einstein/relativistic hydrodynamics equations on a supercomputer with a few simple 
commands. Users need no prior knowledge about either the details of the architecture of a 
supercomputer nor its particular software configuration. The same exact set of SimFactory 
commands can be used on all other supported supercomputers to run the same simulation 
there. 

The Simulation Factory supports and simplifies three kinds of operations: 

1. Remote Access. The actual access commands and authentication methods differ 
between systems, as do the user names that a person has on different systems. Some 
systems are not directly accessible, and one must log in to a particular "trampoline" 
server first. The Simulation Factory hides this complexity. 
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2. Configuring and Building. Building Cactus requires certain software on the system, 

such as compilers, libraries, and build tools. Many systems offer different versions of 
these, which may be installed in non-default locations. Finding a working combination 
that results in efficient code is extremely tedious and requires low-level system 
experience. The Simulation Factory provides a machine database that enables users 
to store and exchange this information. In many cases, this allows people to begin to 
use a new machine in a very short time with just a few, simple commands. 

3. Submitting and Managing Simulations. Many simulations run for days or weeks, 

requiring frequent checkpointing and job re-submission because of short queue run- 
time limits. Simple user errors in these menial tasks can potentially destroy weeks 
of information. The Simulation Factory offers safe commands that encapsulate best 
practices that prevent many common errors and leave a log trail. 

The above features make running simulations on supercomputers much safer and simpler. 
4.4- Kranc 

Kranc [TT6] - lll8] is a Mathematica application which converts a high-level continuum 
description of a PDE into a highly optimized module for Cactus, suitable for running on 
anything from a laptop to the world's largest HPC systems. Many codes contain a large 
amount of complexity, including expanded tensorial expressions, numerical methods, and 
the large amount of "glue" code needed for interfacing a modern HPC application with the 
underlying framework. Kranc absorbs this complexity, allowing the scientist to concentrate 
on writing only the Kranc script which describes the continuum equations. 

This approach brings with it many advantages. With these complicated elements 
factored out, a scientist can write many different Kranc codes, all taking advantage of the 
features of Kranc and avoiding unnecessary or trivial but painstaking duplication. The codes 
might be variations on a theme, perhaps versions which use different sets of variables or 
formulations of the equations, or they could represent completely different physical systems. 
The use of a symbolic algebra package, Mathematica, enables high-level optimizations which 
are not performed by the compiler to be implemented in Kranc. 

Any enhancements to Kranc can be automatically applied to all codes which are 
generated using Kranc. Existing codes have easily benefited from the following features 
added to Kranc after the codes themselves were written: (i) OpenMP parallelization support, 
necessary for efficient use of modern multi-core processors; (ii) support for multipatch 
domains with the Llama |119j code; (iii) automatic generation of vectorized code, where 
the equations are evaluated simultaneously by the processor for two grid points at the same 
time; and (iv) common sub-expression elimination, and various other optimization strategies. 

Within the Einstein Toolkit, the Einstein evolution thorn McLachlan, as well as the 
wave extraction thorn WeylScal4, are both generated using Kranc, and hence support all 
the above features. 
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5. Components 



The Einstein Toolkit uses the modular Cactus framework as its underlying infrastructure. 
A simulation within Cactus could just use one module, but in practice simulations are often 
composed from hundreds of components. Some of these modules provide common definitions 



and conventions (see section 5.1). Others provide initial data (see section 5.2), which may be 



evolved using the different evolution methods for vacuum and matter configurations described 



in sections 5.3 and 5.4, respectively. The thermodynamic properties of fluids are encoded 



in equations of state (see section 5.5). Finally, additional quantities which are not directly 
evolved are often interesting for a detailed analysis of the simulation's results. Modules 



providing commonly used analysis methods are described in section 5.6 



5.1. Base Modules 

Modular designs have proven to be essential for distributed development of complex software 
systems and require the use of well-defined interfaces. Low-level interoperability within the 
Einstein Toolkit is provided by the Cactus infrastructure. One example of this is the usage of 
one module from within another, e.g., by calling a function within another thorn independent 
of programming language used for both the calling and called function. Solutions for technical 
challenges like this can be and are provided by the underlying framework, in this case Cactus. 

However, certain other standards are very hard or impossible to enforce on a technical 
level. Examples for these include the exact definitions of physical variables, their units, and, 
on a more technical level, the variable names used for the physical quantities. Even distinct 
simulation codes typically use very similar scheduling schemes, so conventions describing 
the behavior of the scheduler can help coordinate the order in which functions in different 
modules are called. 

The Einstein Toolkit provides modules whose sole purpose is to declare commonly used 
variables and define their meaning and units. These conditions are not strictly enforced, but 
instead documented for the convenience of the user community. Three of these base modules, 
ADMBase, HydroBase, and TmunuBase, are described in more detail below. 

In the following, we assume that the reader is familiar with the basics of numerical 
relativity and GR hydrodynamics, including the underlying differential geometry and 
tensor analysis. Detailed introductions to numerical relativity have recently been given by 
Alcubierre |120j . Baumgarte & Shapiro |121j . and Centrella et al. |122j . GR hydrodynamics 
has been reviewed by Font [3D] . We set G = c = 1 throughout this paper, and Mq = 1 where 
appropriate. 



5.1.1. ADMBase The Einstein Toolkit provides code to evolve the Einstein equations 

G^"' = 87iT^'', (1) 

where G'^^ is the Einstein tensor, describing the curvature of 4-dimensional spacetime, and 
T'^'^ is the stress-energy tensor. Relativistic spacetime evolution methods used within the 
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Cactus framework employ different formalisms to accomplish this goal, but essentially all 
are based on the 3 + 1 ADM construction [123j, which makes it the natural choice of a 
common foundation for exchange data between modules using different formalisms. In the 
3 + 1 approach, 4-dimensional spacetime is foliated into sequences of spacelike 3-dimensional 
hypersurfaces (slices) connected by timelike normal vectors. The 3 + 1 split introduces 4 
gauge degrees of freedom: the lapse function a that describes the advance of proper time 
with respect to coordinate time for a normal observer§ and the shift vector that describes 
how spatial coordinates change from one slice to the next. 

According to the ADM formulation, the spacetime metric is assumed to take the form 

rfs2 ^ g^^dx^'dx" = (-a^ + P,P')dt^ + 2P,dt dx' + 'jijdx'dx^, (2) 

where gfj_,y, a, f3\ and 'jij are the spacetime 4-metric, lapse function, shift vector, and 
spatial 3-metric, respectively, and we follow the standard relativistic convention where Latin 
letters are used to index 3-dimensional spatial quantities and Greek letters to index 4- 
dimensional spacetime quantities, with the index running from to 3. The remaining 
dynamical component of the spacetime is contained in the definition of the extrinsic curvature 
Kij, which is defined in terms of the time derivative of the metric after incorporating a Lie 
derivative with respect to the shift vector: 

Kij = -^{dt- Cp)-fij. (3) 

The three-metric, extrinsic curvature, lapse function, and shift vector are all declared as 
variables in the ADMBase module, the latter two together with their first time derivatives. 
The variables provided by ADMBase are: 

• The 3-metric tensor, 7^^-: gxx, gxy, gxz,gyy, gyz, gzz 

• The extrinsic curvature tensor, Kij: kxx, kxy, kxz, kyy, kyz, kzz 

• The lapse function, a: alp 

• The shift vector /?*: betax, betay, betaz 

This base module also defines common parameters to manage interaction between 
different modules. Examples are the type of requested initial data or the used evolution 
method. 

The type of initial data chosen for a simulation is specified by the parameters 
initial_data (3-metric and extrinsic curvature), initial_lapse, initial_shif t. The 
time derivatives of the gauge variables (the lapse and shift) are set by the parameters 
initial-dt lapse and initial_dtshif t, respectively. By default, ADMBase initializes the 
3-metric and extrinsic curvature to Minkowski (i.e., 'jij = Sij, the Kronecker delta, and 
Kij = 0), the shift to zero, and the lapse to unity. Initial data thorns override these defaults 
by extending the parameters. 

Analogous to specifying initial data, evolution methods are chosen by the 
parameters evolutionjnethod (3-metric and extrinsic curvature), lapse_evolutionjnethod, 

§ A normal observer follows a worldline tangent to the unit normal on the 3-hypersurface. 
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shif t_evolution_method, dtlapse_evolution_method and dtshif t.evolutionjiiethod. 
ADMBase does not evolve the 3-metric or extrinsic curvature, and holds the lapse and shift 
static. Evolution thorns extend the ranges of these parameters and contain the evolution 
code. 

The variables defined in ADMBase typically are not used for the actual evolution of 
the curvature. Instead, it is expected that every evolution module converts its internal 
representation to the form defined in ADMBase after each evolution step. This procedure 
enables modules which perform analysis on the spacetime variables to use the ADMBase 
variables without direct dependency on any of the existing curvature evolution methods. 

5.1.2. HydroBase Similar to ADMBase, the module HydroBase defines a common basis 
for interactions between modules of a given evolution problem, in this case relativistic 
hydrodynamics. HydroBase extends the Einstein Toolkit to include an interface within which 
magnetohydrodynamics may work. HydroBase's main function is to store variables which are 
common to most if not all hydrodynamics codes solving the Euler equations, the so-called 
primitive variables. These are also the variables which are needed to couple to a spacetime 
solver, and often by analysis thorns as well. As with ADMBase, the usage of a common set 
of variables by different hydrodynamics codes creates the possibility of sharing parts of the 
code, e.g., initial data solvers or analysis routines. HydroBase also defines commonly needed 
parameters and schedule groups for the main functions of a hydrodynamics code. 

HydroBase uses a set of conventions known as the Valencia formulation |124H126] . In 
particular, HydroBase defines the primitive variables (see [lO] for details): 

• rho: rest mass density p 

• press: pressure P 

• eps: internal energy density e 

• vel [3] : contravariant fluid three velocity f ' defined as 







+ 



(4) 



aw 



a 



in terms of the four- velocity m^, lapse, and shift vector . 

Y_e: electron fraction Ye 

temperature: temperature T 

entropy: specific entropy per particle s 

Bvec [3] : contravariant magnetic field vector defined as 




(5) 



in terms of the dual F*^'^ = ^e^'^°'^Fai3 to the Faraday tensor and the unit normal of the 
foliation of spacetime ra^ = a~^[l, —P^]^- 
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HydroBase also sets up scheduling blocks that organize the main functions which modules 

of a hydrodynamics code may need. All of those scheduling blocks arc optional, but when 
used they simplify existing codes and make them more interoperable. HydroBase itself does 
not schedule routines inside most of the groups that it provides. Currently the scheduling 

blocks are: 

• Initializing the primitive variables 

• Converting primitive variables to conservative variables 

• Calculating the right hand side (RHS) in the method of hues (MoL) 

• Setting and updating an excision mask 

• Applying boundary conditions 

Through these, the initiation of the primitive variables, methods to recover the 
conservative variables, and basic atmosphere handling can be implemented in different thorns 
while allowing a central access point for analysis thorns. 

5.1.3. TmunuBase In the Einstein Toolkit, we typically choose the stress energy tensor T^^ 
to be that of an ideal relativistic fluid, 

T'*'^ = phui'u'' - g^^'P , (6) 

where p, u^, and g^" are defined above, and h — 1+e+P/pis the relativistic specific enthalpy. 

The thorn TmunuBase provides grid functions for the stress-energy tensor T^,^ as well as 
schedule groups to manage when T^^ is calculated. In a simulation, many different thorns may 
contribute to the stress-energy tensor and this thorn allows them to do so without explicit 
interdependence. The resulting stress-energy tensor can then be used by the spacetime 
evolution thorn (again without explicit dependence on any matter thorns). When thorn MoL 
is used for time evolution this provides a high-order spacetime-matter coupling. 

The grid functions provided by TmunuBase are: 

• The time component Tqo: eTtt 

• The mixed components Toj: eTtx, eTty, eTtz 

• The spatial components Tj^: eTxx, eTxy, eTxz, eTyy, eTyz, eTzz 

In addition, the grid scalar stress_energy_state has the value 1 if storage is provided for 
the stress-energy tensor and if not. 

Thorn ADMCoupling provides a similar (but older) interface between space-time and 
matter thorns. However, since it is based on an include file mechanism it is more complicated 
to use. We recommend all new matter thorns to use TmunuBase instead. 

5.2. Initial Data 

The Einstein Toolkit contains many modules used to generate initial data for GR simulations, 
including both vacuum and hydrodynamic configurations. These include modules used 
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primarily for testing of various components, as well as physically motivated configurations 
that describe, for example, single or binary BHs and/or NSs. Many of the modules are 
self-contained, consisting of either all the code to generate exact initial solutions or the 
numerical tools required to construct solutions known semi-analytically. Others, though, 
require the installation of numerical software packages that are included in the toolkit as 
external libraries. One example is the TwoPunctures module |127] — commonly used in 
numerical relativity to generate BH-BH binary initial data — which makes use of the GNU 
Scientific Library [GSL; |128l 1129] ]. Several modules have also been implemented to read in 
data files generated by the Lorene code \130\ I131j . 

Initial data setup is in most cases clearly separated from the evolution that follows. 
Typically, initial data routines provide the data in terms of the quantities defined in the 
Base modules (see section 5.1), and the evolution modules will later convert these quantities 
to forms used for the evolution. For example, an initial data module must supply Qij, the 
spatial 3-metric, and Kij, the extrinsic curvature. The conversion between the physical 
metric and extrinsic curvature and conformal versions of these is handled solely within 
evolution modules, which are responsible for calculating the conformally related three metric 
7jj, the conformal factor ip, the conformal traceless extrinsic curvature Aij and the trace 
of the extrinsic curvature K, as well as initializing the BSSN variable F* should that be 
the evolution formalism chosen (see section 5.3 for definitions of these). Optionally, many 
initial data modules also supply values for the lapse and shift vector and in some cases their 
time derivatives. It is important to note, though, that many dynamical calculations run 
better from initial gauge choices set by ansatz rather than those derived from stationary 
approximations that are incompatible with the gauge evolution equations. In particular, 
conformal thin-sandwich initial data for binaries include solutions for the lapse and shift that 
are frequently replaced by simpler analytical values that lead to more circular orbits under 
standard "moving puncture" gauge conditions (see, e.g., fI2 \ I132j and other works). 

We turn our attention next to a brief discussion of the capabilities of the aforementioned 
modules as well as their implementation. 



5.2.1. Simple Vacuum initial data Vacuum spacetime tests in which the constraint 
equations are explicitly violated are provided by IDConstraintViolate and Exact, a set 
of exact spacetimes in various coordinates including Lorentz-boosted versions of traditional 
solutions. Vacuum gravitational wave configurations can be obtained by using either 
IDBrillData, providing a Brill wave spacetime |133 j: or IDLinearWaves, for a spacetime 
containing a linear gravitational wave. Single BH configurations include IDAnalyticBH 
which generates various analytically known BH configurations; as well as IDAxibrillBH, 
IDAxiOddBrillBH, DistortedBHIVP and RotatingDBHIVP, which introduce perturbations to 
isolated BHs. 



5.2.2. Hydrodynamics Tests Initial data to test different parts of hydrodynamics evolution 
systems are provided by GRHydro_InitData. This module includes several shock tube 
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problems that may be evolved in any of the Cartesian directions or diagonally. All of 
these have been widely used throughout the field to evaluate a diverse set of solvers |134j . 
Conservative to primitive variable conversion and vice versa are also supported, as are tests 



to check on the reconstruction of hydrodynamical variables at cell faces (see Sec. 5.4 for more 
on this). Along similar lines, the Hydro_InitExcision module sets up masks for different 
kinds of excised regions, which is convenient for hydrodynamics tests. 

5.2.3. TwoPunctures: Binary Black Holes and extensions A substantial fraction of the 
published work on the components of the Einstein toolkit involves the evolution of BH- 
BH binary systems. The most widely used routine to generate initial data for these is the 
TwoPunctures code (described originally in |127] ) which solves the binary puncture equations 
for a pair of BHs |135j . To do so, one assumes the extrinsic curvature for each BH corresponds 
to the Bowen-York form |136j : 

+ ^.{e'^'st^NiN^ + E^''sjr'^NiN'), (7) 

where the sub/superscript (m) refers to the contribution from BH m = 1, 2; the 3-momentum 
is p*; the BH spin angular momentum is Sf, the conformal 3-metric 7*-' is assumed to be flat, 
i.e. = rjij, and A^* = x^/r is the Cartesian normal vector relative to the position of each BH 
in turn. This solution automatically satisfies the momentum constraint, and the Hamiltonian 
constraint may be written as an elliptic equation for the conformal factor, defined by the 
condition gij = ip'^'jij or equivalently ip = (det \gij\Y^^'^- 

+ Ik'^K,^^-' = (8) 

o 

Decomposing the conformal factor into a singular analytical term and a regular term u, such 
that 

mi 772,2 1 / N 

2ri 2r2 * 

where mi, m2 and ri, r2 are the mass of and distance to each BH, respectively, and \E' is 
defined by the equation itself, the Hamiltonian constraint may be written as 
'1 



Au + 



I ''J 



;i + ^n)-^ (10) 



subject to the boundary condition m — )■ 1 as r — )■ 00. In Cartesian coordinates, the function 
u is infinitely different iable everywhere except the two puncture locations. TwoPunctures 
resolves this problem by performing a coordinate transformation modeled on confocal 
elliptical/hyperbolic coordinates. This transforms the spatial domain into a finite cube 
with the puncture locations mapped to two parallel edges, as can be seen in figure [3] The 
coordinate transformation is: 

_ /12 + 1 2B 
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Figure 3. Example of a TwoPunctures coordinate system for BH-NS binary initial data 



, 2A 1-52 

' = \-A^l + B^ ''^'^ (11) 
which maps TZ^ onto < A < 1 (the elhptical quasi-radial coordinate), — 1 < -B < 1 (the 
hyperbohc quasi-latitudinal coordinate), and < < 27r (the longitudinal angle). Since 
u is smooth everywhere in the interior of the remapped domain, expansions into modes in 
these coordinates are spectrally convergent and thus capable of extremely high accuracy. In 
practice, the field is expanded into Chebyshev modes in the quasi-radial and quasi-latitudinal 
coordinates, and into Fourier modes around the axis connecting the two BHs. The elliptic 
solver uses a stabilized biconjugate gradient method to achieve rapid solutions and to avoid 
ill-conditioning of the spectral matrix. 



5.2.4- Lorene-based binary data The ET contains three routines that can read in publicly 
available data generated by the Lorene code |13Ul 1131] . though it does not currently 
include the capability of generating such data from scratch. For a number of reasons, such 
functionality is not truly required; in particular, Lorene is a serial code and to call it as an 
ET initial data generator saves no time. Also, it is not guaranteed to be convergent for an 
arbitrary set of parameters; thus the initial data routine itself may never finish its iterative 
steps. Instead, recommended practice is to let Lorene output data into files, and then read 
those into ET at the beginning of a run. 

Lorene uses a multigrid spectral approach to solve the conformal thin-sandwich equations 
for binary initial configurations |132j and a single-grid spectral method for rotating stars. 
For binaries, five elliptic equations for the shift, lapse, and conformal factor are written 
down and the source terms are divided into pieces that are attributed to each of the two 
objects. Matter source terms are ideal for this split, since they are compactly supported, while 
extrinsic curvature source terms are spatially extended but with sufficiently rapid fallofT at 
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Figure 4. Example of a Lorene multi-domain coordinate system for binary initial data. The 
outermost, compactified domain extending to spatial infinity is not shown. 



large radii to yield convergent solutions. Around each object, a set of nested spheroidal sub- 
domains (see figure |4]) is constructed to extending through all of space, with the outermost 
domain incorporating a compactification to allow it to extend to spatial infinity. Within 
each of the nested sub-domains, fields are decomposed into Chebyshev modes radially and 
into spherical harmonics in the angular directions, with elliptic equation solving reduced to 
a matrix problem. The nested sub-domains need not be perfectly spherical, and indeed one 
may modify the outer boundaries of each to cover any convex shape. For NSs, this allows 
one to map the surface of a particular sub-domain to the NS surface, minimizing Gibbs error 
there. For BHs, excision boundary conditions are imposed at the horizon. To read a field 
solution describing a given quantity onto a Cactus-based grid, one must incorporate the data 
from each star's domains at that point. 

Meudon_Bin_BH can read in BH-BH binary initial data described in [137j . while 
Meudon_Bin_NS handles binary NS data from jl31j . Meudon_Mag_NS may be used to read 
in magnetized isolated NS data |13U] . 

5.2.5. TOVSolver The TOVSolver routine in the ET solves the standard TOV equations 
|138[ 1139] expressed using the Schwarzschild (areal) radius r in the interior of a spherically 
symmetric star in hydrostatic equilibrium: 



where e = p(l + e) is the energy density of the fiuid, including the internal energy 
contributionf, m is the gravitational mass inside a sphere of radius r, and $ the logarithm 




dm 
dr 

dr 




(12) 



f We note that since different application thorns may define their own local variables, the energy density 
is referred to as rho within TOVSolver, as the projected energy density E, defined in Sec. |5.3[ is within 
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of the lapse. The routine also supplies the analytically known solution in the exterior, 
P = TOV_atmo sphere, 

m = M, 

$ = ^log(l - 2M/r) (13) 

where TOV_atmo sphere is a parameter used to define the density of the ambient atmosphere. 
Since the isotropic radius f is the more commonly preferred choice to initiate dynamical 
calculations, the code then transforms all variables into isotropic coordinates, integrating the 
radius conversion formula 

d{\og{f/r)) r^/^ — (r — 2m)-^/^ 



(14) 



dr r{r — 2m) 

subject to the boundary condition that in the exterior, 

f = ^ (^Vr2 - 2Mr + r - M j 

1 + ^) . (15) 

handling with some care the potentially singular terms that appear at the origin. 

To facilitate the construction of stars in more complicated dynamical configurations, one 
may also apply a uniform velocity to the NS, though this does not affect the ODE solution 
nor the resulting density profile. 

5.3. Spacetime Curvature Evolution 

The Einstein Toolkit curvature evolution code McLachlan |112l I140j is auto-generated from 



tensor equations via Kranc (Sec. 4.4). It implements the Einstein equations in a 3 + 1 split 
as a Cauchy initial boundary value problem |141] . For this, the Baumgarte-Shapiro-Shibata- 
Nakamura (BSSN) conformal-tracefree reformulation [9], [TOl I142j of the original Arnowitt- 
Deser-Misner (ADM) formalism |123] is employed. McLachlan uses fourth-order accurate 
finite differencing for the spacetime variables and adds a fifth-order Kreiss-Oliger dissipation 
term to remove high frequency noise. The evolved variables are the conformal factor $, 
the conformal 3-metric 7jj, the trace K of the extrinsic curvature, the trace free extrinsic 
curvature Aij and the conformal connection functions F*. These are defined in terms of the 
standard ADM 4- metric gij, 3-metric jij, and extrinsic curvature Kij by: 

1 



log 



^2 



(16) 



7,, ^e-^S,,, (17) 

K =g'^K,j, (18) 

McLachlan and a few other thorns. Similar ambiguities exist for other commonly used variable names, 
particularly (f> and a. 
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1 



-A(j> 

The evolution equations are then: 



doK 

do/3' 
doB' 

do(p 



— — e 



-40 



+ a I A'^Aij + -K 



aS 



^a^G{a,(j),x'')B' 
6 6 

I. 2 . 



doAij 



aRij + aRfj — DiDja + 49(j0 • i^j)Q; 



(19) 
(20) 

(21) 

(22) 

(23) 
(24) 

(25) 
(26) 

(27) 



- 2A'W^a + 2a 



O771 _ 

(m - l)9fcA^^ - -^D'K 



+ m(f ^^i'^' + Q>A'^dj(t)) 



-S\ 



(28) 



with the momentum constraint damping constant set to m = 1. The stress energy tensor 
Tj^u is incorporated via the projections 



E = — (Too - + I3'P^T'^) 

S = 9"% 



Si = -- (To, - 13%) 



(29) 
(30) 

(31) 



We have introduced the notation = dt — ^^dy AH quantities with a tilde involve the 
conformal 3-metric 7jj, which is used to raise and lower indices. In particular, Di and 
rfj- refer to the covariant derivative and the Christoffel symbols with respect to 7jj. The 
expression [• • denotes the trace-free part of the expression inside the parentheses, and 
we define the Ricci tensor contributions as: 

= - lf%daij + %i^^J)r' - Timdn"' + 1 (2f '^(.f + f '^.f (32) 

Rt = - 2 A - '^lijD''Dk(t> + Dj(t> - ^^i^D^ Dk(t>- (33) 



This is a so-called 0- variant of BSSN. The evolved gauge variables are lapse ct, shift and 
a quantity B^ related to the time derivative of the shift. The gauge parameters /, G, H, and 
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7] are determined by our choice of a 1 + log |143] slicing: 

f{a,^,x>^)=2/a (34) 

Koix^") = (35) 
and F-driver shift condition |143] : 

G{a,(f),x^) = (3/4) (36) 

H{a,(l),x>') = exp{40} (37) 

r]{B\ a, x^") = (1/2) B'q{r). (38) 

The expression g(r) attenuates the F-driver depending on the radius as described below. 

The F-driver shift condition is symmetry-seeking, driving the shift /3* to a state that 
renders the conformal connection functions F* stationary. Of course, such a stationary state 
cannot be achieved while the metric is evolving, but in a stationary spacetime the time 
evolution of the shift /3* and thus that of the spatial coordinates a;* will be exponentially 



damped. This damping time scale is set by the gauge parameter rj (see (38)) which has 
dimension 1/T (inverse time). As described in |144l I145j . this time scale may need to be 
adapted in different regions of the domain to avoid spurious high-frequency behavior in 
regions that otherwise evolve only very slowly, e.g., far away from the source. 

Here we use the simple damping mechanism described in (12) of |145] . which is defined 



as: 



_ . 1 for r < i? (near the origin) 
' R/r for r > R (far away) 

with a constant R defining the transition radius between the interior, where g ~ 1, and the 
exterior, where q falls off as 1/r. A description of how q appears in the gauge parameters 



may be found in (38). 



5.3.1. Initial Conditions Initial conditions for the ADM variables Qij, Kij, lapse a, and 



shift /3* are provided by the initial data routines discussed in Sec. |5.2[ From these the BSSN 
quantities are calculated via their definitions, setting = 0, and using cubic extrapolation 
for F* at the outer boundary. This extrapolation is necessary since the F* are calculated from 
derivatives of the metric, and one cannot use centered finite differencing stencils near the 
outer boundary. 

The extrapolation stencils distinguish between points on the faces, edges, and corners of 
the grid. Points on the faces are extrapolated via stencils perpendicular to that face, while 
points on the edges and corners are extrapolated with stencils aligned with the average of the 
normals of the adjoining faces. For example, points on the {+x, +y) edge are extrapolated 
in the (1,1,0) direction, while points in the {+x,+y + z) corner are extrapolated in the 
(1, 1, 1) direction. Since several layers of boundary points have to be filled for higher order 
schemes (e.g., three layers for a fourth order scheme), one proceeds outwards starting from 
the innermost layer. Each subsequent layer is then defined via the points in the interior and 
the previously calculated layers. 
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5.3.2. Boundary Conditions During time evolution, a Sommerfeld-type radiative boundary 
condition is applied to all components of the evolved BSSN variables as described in |142] . 
The main feature of this boundary condition is that it assumes approximate spherical 
symmetry of the solution, while applying the actual boundary condition on the boundary of 
a cubic grid where the face normals are not aligned with the radial direction. This boundary 
condition defines the right hand side of the BSSN state vector on the outer boundary, which 
is then integrated in time as well so that the boundary and interior are calculated with the 
same order of accuracy. 

The main part of the boundary condition assumes that one has an outgoing radial wave 
with some speed Vq: 

X = X„+^±^, (40) 

r 

where X is any of the tensor components of evolved variables, Xq the value at infinity, and 
u a spherically symmetric perturbation. Both Xq and Vq depend on the particular variable 
and have to be specified. This implies the following differential equation: 

dtX = - v'd,X - vo , (41) 

r 

where f* = vox^/r. The spatial derivatives di are evaluated using centered finite differencing 
where possible, and one-sided finite differencing elsewhere. Second order stencils are used in 
the current implementation. 

In addition to this main part, it is also necessary to account for those parts of the 
solution that do not behave as a pure wave, e.g.. Coulomb type terms caused by infall of 
the coordinate lines. It is assumed that these parts decay with a certain power p of the 
radius. This is implemented by considering the radial derivative of the source term above, 
and extrapolating according to this power-law decay. 

Given a source term {dtX), one defines the corrected source term (dtX)* via 

{dtXy = (dtX) + ( n^d.idtX) , (42) 

where n* is the normal vector of the corresponding boundary face. The spatial derivatives di 
are evaluated by comparing neighboring grid points, corresponding to a second-order stencil 
evaluated in the middle between the two neighboring grid points. Second-order decay is 
assumed, hence p = 2. 

As with the initial conditions above, this boundary condition is evaluated on several 
layers of grid points, starting from the innermost layer. Both the extrapolation and radiative 
boundary condition algorithms are implemented in the publicly available NewRad component 
of the Einstein Toolkit. 

This boundary condition is only a coarse approximation of the actual decay behavior 
of the BSSN state vector, and it does not capture the correct behavior of the evolved 
variables. However, one finds that this boundary condition leads to stable evolutions if 
applied sufficiently far from the source. Errors introduced at the boundary (both errors 
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in the geometry and constraint violations) propagate inwards with the speed of hght |140] . 
Gauge changes introduced by the boundary condition, which are physically not observable, 
propagate faster, with a speed up to \/2 for the gauge conditions used in McLachlan. 

5.4- Hydrodynamics Evolution 

Hydrodynamic evolution in the Einstein Toolkit is designed so that it interacts with the metric 
curvature evolution through a small set of variables, allowing for maximum modularity in 
implementing, editing, or replacing either evolution scheme. 

The primary hydrodynamics evolution routine in the Einstein Toolkit is GRHydro, a code 
derived from the public Whisky code [OU I146H148] designed primarily by researchers at AEI 
and their collaborators. It includes a high resolution shock capturing (HRSC) scheme to 
evolve hydrodynamic quantities, with several different reconstruction methods and Riemann 
solvers, as we discuss below. In such a scheme, we define a set of "conserved" hydrodynamic 
variables, defined in terms of the "primitive" physical variables such as mass and internal 
energy density, pressure, and velocity. Wherever derivatives of hydrodynamic terms appear in 
the evolution equations for the conserved variables, they are restricted to appear only inside 
divergence terms (referred to as fluxes) and never in the source terms. By calculating fluxes at 
cell faces, we may obtain a consistent description of the inter-cell values using reconstruction 
techniques that account for the fact that hydrodynamic variables are not smooth and may 
not be finite differenced accurately. All other source terms in the evolution equations may 
contain only the hydrodynamic variables themselves and the metric variables and derivatives 
of the latter, since the metric must formally be smooth and thus differentiable using finite 
differencing techniques. Summarizing these methods briefly, the following stages occur every 
timestep: 

• The primitive variables are "reconstructed" at cell faces using shock-capturing 
techniques, with total variation diminishing (TVD), piecewise parabolic (PPM), and 
essentially non-oscillatory (ENO) methods currently implemented. 

• A Riemann problem is solved at each cell face using an approximate solver. Currently 
implemented versions include HLLE (Harten-Lax-van Leer-Einfeldt), Roe, and Marquina 
solvers. 

• The conserved variables are advanced one timestep, and used to recalculate the new 
values of the primitive variables. 

We discuss the GRHD formalism, the stages within a timestep, and the other aspects of the 
code below, noting that the documentation included in the released version is quite extensive 
and covers many of these topics in substantially more detail. 

5.4-i- Ideal general relativistic hydrodynamics (GRHD) The equations of ideal GR 
hydrodynamics evolved by GRHydro are derived from the local GR conservation laws of mass 
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and energy-momentum: 

V^J^ = 0, V^T'^" = , (43) 

where denotes the covariant derivative with respect to the 4-metric, and = p-u^ is the 
mass current. 

The 3-velocity f * may be calculated in the form 

where W = (1 — v^Vi) is the Lorentz factor. The contravariant 4- velocity is then given 
by: 



u" = —, u' = W [v'- — ] , (45) 
a \ a J 

and the covariant 4-velocity is: 

no = W{v'f3i -a), = Wvi . (46) 

The GRHydro evolution scheme is a first-order hyperbolic flux-conservative system for 
the conserved variables D, S**, and r, which may be defined in terms of the primitive variables 
p, e, v\ such that: 

D = ^pW, (47) 
= ^phW^v\ (48) 
T = ^{phW^ ~ P) - D , (49) 

where 7 is the determinant of 7jj. The evolution system then becomes 



with 



= a [Dv\ Sjv' + 5;P, Tv' + Pv'], 



S = a 



a 



Q^lna ^ rpfiv-pQ 



(51) 



Here, n* = f * — /3Yq; and F^^ are the 4-Christoffel symbols. The time integration and 
coupling with curvature are carried out with the Method of Lines |149] . The expressions for 
S are calculated in GRHydro by using the definition of the extrinsic curvature to avoid any 
time derivatives whatsoever, as discussed in detail in the code's documentation, following a 
suggestion by Mark Miller based on experience with the GR3D code. 
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5.4-2. Reconstruction techniques In order to calculate fluxes at cell faces, we first must 
calculate function values on either side of the face. In practice, reconstructing the primitive 
variables yields more stable and accurate evolutions than reconstructing the conservatives. 
In what follows, we assume a Cartesian grid and loop over faces along each direction in 
turn. We define q!['^i/2 to be the value of a quantity q on the left side of the face between 
qi = q{xi, y, z) and g^+i = y, z), where Xi is the ith point in the x-direction, and (?f^i/2 

the right side of the same face. 

For total variation diminishing (TVD) methods, we let: 



?i+i/2 — Q^i -1 2 ' ^*+i/2 ~ 2 ^ ' 



where f{qi) is a slope-limited gradient function, typically determined by the values of g^+i — gj 
and qi — ^j-i, with a variety of different forms of the slope limiter available. In practice, all 
try to accomplish the same task of preserving monotonicity and removing the possibility of 
spuriously creating local extrema. Implemented methods include minmod, superbee |150] . 
and monotonized central |151j . 

The piecewise parabolic method (PPM) is a multi-step method based around a quadratic 
fit to nearby points interpolated to cell faces |152] , for which q^ and q^ are generally equivalent 
except near shocks and local extrema. The version implemented in GRHydro includes the 
steepening and flattening routines described in the original PPM papers, with a simplified 
dissipation procedure. Essentially non-oscillatory (ENO) methods use a divided differences 
approach to achieve third-order accuracy via polynomial interpolation |153[ 1154] . 

Both ENO and PPM yield third-order accuracy for smooth monotonic functions, whereas 
TVD methods typically yield second-order accurate values. Regardless of the reconstruction 
scheme chosen, all of these methods reduce to first order near local extrema and shocks. 



5.4-3. Riemann solvers The Riemann problem involves the solution of the equation 

dtq + d,f\q) = Q (53) 

at some point X representing a discontinuity between constant states. The exact solution can 
be quite complicated, involving five different waves with different characteristic speeds for a 
hydrodynamic problem (8 for GRMHD), so GRHydro implements three different approximate 
solvers to promote computational efficiency. In each case, the solution takes a self-similar 
form q{C}-, where ^ = (x — X)/t represents the characteristic speed from the original shock 
location to the point in question in space and time. 

The simplest method implemented is the Harten-Lax-van Leer-Einfeldt solver |155[ I156j 
(HLL or HLLE, depending on the reference), which uses a two wave approximation to 
calculate the evolution along the shock front. With ^_ and the most negative and most 
positive wave speeds present on either side of the interface, the solution g(^) is assumed to 
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take the form 

' if e < 



q={ if < e < e+ (54) 
if e>^+, 

with the intermediate state g^, given by 

g* — £• £• ■ y^^) 

The numerical flux along the interface takes the form 



where 



= min(0,^), e+ = max(0,e+). (57) 



It is these flux terms that are then used to evolve the hydrodynamic quantities. 

The Roe solver |157] involves linearizing the evolution system for the hydrodynamic 



evolution, defining the Jacobian matrix A = ^ (see (53)), and working out the eigenvalues 



A* and left and right eigenvectors, k and , assumed to be orthonormalized so that k-r^ = 5^. 
Defining the characteristic variables Wi = It ■ q, the characteristic equation becomes 

dtw + Ad^w = (58) 

with A the diagonal matrix of eigenvalues. Letting Awi = — wf = U ■ {q^ — q^) represent 
the differences in the characteristic variables across the interface, the Roe flux is calculated 

as 

f{q) = I (/(g^) + fiq"") - E IM^^^r^) (59) 

where the eigenvector appearing in the summed term are evaluated for the approximate Roe 
average flux g^oe = |(g^ + g^). The Marquina flux routines use a similar approach to the Roe 
method, but provide a more accurate treatment for supersonic flows (i.e., those for which the 
characteristic wave with ^ = is within a rarefaction zone) [158[ 1159] . 



5. 4-4- Conservative to primitive conversion In order to invert eqs. (47) - (49 ), solving for the 
primitive variables based on the values of the conservative ones, GRHydro uses a 1-dimensional 
Newton- Raphson approach that solves for a consistent value of the pressure. Defining the 
(known) undensitized conservative variables D = -D/a/7 = pW, = S"" / ^ = phW^v^ and 
f = = phW^ — P — D, as well as the auxiliary quantities Q = phW'^ = f + D + P and 

= 'jijS^S^ = {phW)'^{W'^ — 1), the former of which depends on P and the latter of which 



is known, we find that y — S"^ = phW and thus 
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W = , ^ (61) 



Q'^-S'^-PW -D 

■ (62) 



D 

Given the new values of p and e, one may then find the residual between the pressure and 
P(p, e) and perform the Newton-Raphson step, so long as the values of ^ and ^ are known. 

5.4-5. Atmospheres, boundaries, and other code details GRHydro uses an atmosphere, or 
extremely-low density floor, to avoid problems involving sound speeds and conservative- 
to-primitive variable conversion near the edges of matter distributions. The floor density 
value may be chosen in either absolute (rho_abs_min) or relative (rho_rel_inin) terms. The 
atmosphere is generally assumed to have a specified polytropic EOS, regardless of the EOS 
describing the rest of the matter within the simulation. Whenever the numerical evolution 
results in a grid cell where conservative to primitive variable conversion yields negative values 
of either p or e, the cell is reassigned to the atmosphere, with zero velocity. 

At present, only flat boundary conditions are supported for hydrodynamic variables, 
since it is generally recommended that the outer boundaries of the simulation be placed 
far enough away so that all cells near the edge of the computational domain represent the 
atmosphere. 

GRHydro has the ability to advect a set of passive scalars, referred to as "tracers", as 
well as the electron fraction of a fluid, under the assumption that each tracer X follows the 
conservation law 

dt{DX) + di{av'DX) = 0. (63) 
5.5. Equations of State 

An equation of state connecting the primitive state variables is needed to close the system 
of GR hydrodynamics equations. The module EOS_Omni provides a unified general equation 
of state (EOS) interface and back-end for simple analytic and complex microphysical EOSs. 
The polytropic EOS 

p = (64) 

where K is the polytropic constant and F is the adiabatic index, is appropriate for adiabatic 
(= isentropic) evolution without shocks. When using the polytropic EOS, one does not need 
to evolve the total fluid energy equation, since the speciflc internal energy e is flxed to 

Note that the adiabatic index F = o?lnP/cilnp is related to the frequently used polytropic 
index n via n = 1/(F — 1). 
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The gamma-law EOS J, 

P = (r - l)pe , (66) 

allows for non-adiabatic flow but still assumes flxed microphysics, which is encapsulated in 
the constant adiabatic index T. This EOS has been used extensively in simulations of NS-NS 
and BH-NS mergers. 

The hybrid EOS, flrst introduced by |160j . is a 2-piecewise polytropic with a thermal 
component designed for the application in simple models of stellar collapse. At densities 
below nuclear density, a polytropic EOS with F = Fi ^ 4/3 is used. To mimic the stiffening 
of the nuclear EOS at nuclear density, the low-density polytrope is fltted to a second polytrope 
with F = F2 ^ 2. To allow for thermal contributions to the pressure due to shock heating, a 
gamma-law with F = F^h is used. The full EOS then reads 

P = L:il^i^pr,-r r _ (F^h - 1)(F - FQ r 

+ (Fth - l)pe . (67) 

In this, the total speciflc internal energy e consists of a polytropic and a thermal contribution. 
In iron core collapse, the pressure below nuclear density is dominated by the pressure of 
relativistically degenerate electrons. For this, one sets K = 4.897 x 10^^ [cgs] in the above. 
The thermal index Fth is usually set to 1.5, corresponding to a mixture of relativistic (F = 4/3) 
and non-relativistic (F = 5/3) gas. Provided appropriate choices of EOS parameters (e.g., 
|161] ). the hybrid EOS leads to qualitatively correct collapse and bounce dynamics in stellar 
collapse. 

EOS.Omni also integrates the nuc_eos driver routine, which was first developed for the 
GRID code ^49] for tabulated microphysical finite-temperature EOS which assume nuclear 
statistical equilibrium (NSE). nuc_eos handles EOS tables in HDF5 format which contain 
entries for thermodynamic variables X = X{p,T,Ye), where T is the matter temperature 
and Yf. is the electron fraction. nuc_eos also supports calls for X = X{p, e. Ye) and carries 
out a Newton iteration to find T{p, e. Ye). For performance reasons, nuc.eos employs simple 
tri-linear interpolation in thermodynamic quantities and thus requires finely spaced tables 
to maintain thermodynamic consistency at an acceptable level. EOS tables in the format 
required by nuc_eos are freely available from http://stellarcollapse.org. 



5.6. Analysis 

It is often beneficial and sometimes necessary to evaluate analysis quantities during the 
simulation rather than post-processing variable output. Beyond extracting physics, these 
quantities are often used as measures of how accurately the simulation is progressing. In the 
following, we describe the common quantities available through Einstein Toolkit modules, and 
how different modules approach these quantities with differing assumptions and algorithms. 
The most common analysis quantities provided fall broadly into several different categories, 

I For historic reasons, this EOS is referred to as the "ideal fluid" EOS in GRHydro. 
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including horizons, masses and momenta, and gravitational waves. Note that several modules 
bridge these categories and some fall outside them, including routines to perform constraint 
monitoring and to compute commonly used derived spacetime quantities. The following 
discussion is meant as an overview of the most common tools rather than an exhaustive list 
of the functionality provided by the Einstein Toolkit. In most cases, the analysis modules 



work on the variables stored in the base modules discussed in Sec. 5J^, ADMBase, TmunuBase, 
and HydroBase, to be as portable as possible. 

5.6.1. Horizons When spacetimes contain a BH, localizing its horizon is necessary for 
describing time-dependent quasi-local measures such as its mass and spin. The Einstein 
Toolkit provides two modules — AHFinder and AHFinderDirect — for locating the apparent 
horizons (AHs) defined locally on a hypersurface. The module EHFinder is also available to 
search an evolved spacetime for the globally defined event horizons. 

EHFinder [162] evolves a null surface backwards in time given an initial guess (e.g., the 
last apparent horizon) which will, in the vicinity of an event horizon, converge exponentially 
to its location. This must be done after a simulation has already evolved the initial data 
forward in time with enough 3D data written out that the full 4-metric can be recovered at 
each timestep. 

In EHFinder, the null surface is represented by a function = which is required 

to satisfy the null condition g°'^dafdpf = 0. In the standard numerical 3+1 form of the 
metric, this null condition can be expanded out into an evolution equation for / as 



dtf = mf - V<^'r'djdjf (68) 

where the roots are chosen to describe outgoing null geodesies. The function / is chosen such 
that it is negative within the initial guess of the horizon and positive outside, initially set to a 
distance measure from the initial surface guess /(to, a;*) = a/ (x* — xl){xi — Xi(o)) — To. There 
is a numerical problem with the steepening of V/ during the evolution, so the function / is 
regularly re-initialized during the evolution to satisfy |V/| ~ 1. This is done by evolving 

^ = A (|V/I -1) (69) 

for some unphysical parameter A until a steady state has been reached. As the isosurface 
/ = converges exponentially to the event horizon, it is useful to evolve two such null surfaces 
which bracket the approximate position of the anticipated event horizon to further narrow 
the region containing the event horizon. 

However, event horizons can only be found after the full spacetime has been evolved. 
It is often useful to know the positions and shapes of any BH on a given hypersurface for 
purposes such as excision, accretion, and local measures of its mass and spin. The Einstein 
Toolkit provides several algorithms of varying speed and accuracy to find marginally trapped 
surfaces, of which the outermost are AHs. All finders make use of the fact that null geodesies 
have vanishing expansion on an AH which, in the usual 3+1 quantities, can be written 

e = V^n* + Ki^n'v? -K = (70) 
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where n* is the unit outgoing normal to the 2-surface. 

The module AHFinder provides two algorithms for locating AHs. The minimization 
algorithm |163] finds the local minimum of — QoY'P'S corresponding to a surface of 

constant expansion Go, with Go = corresponding to the AH. For time-symmetric data, the 
option exists to find instead the minimum of the surface area, which in this case corresponds 
to an AH. An alternative algorithm provided by AHFinder, the fiow algorithm |164] . on 
which EHFinder is also based. Defining a surface as a level set /(x*) = r — h{6,(j)) = 0, 
and introducing an unphysical timelike parameter A to parametrize the fiow of h towards a 



solution, ( 70 ) can be rewritten 



' max V 

2 • 



where p is a strictly positive weight, L is the Laplacian of the 2D metric, and a, (3, and £max 
are free parameters. Decomposing h{6, (p) onto a basis of spherical harmonics, the coefficients 
agm evolve iteratively towards a solution as 

(n+l) _ (n) Q + /3£max (^max + 1) ( nP\\^"'^ (70\ 

" '"^ W(Cax + l)(l + W+l)/«)^^''^^- ^''^ 

The AHFinderDirect module |165j is a faster alternative to AHFinder. Its approach is 



to view (70) as an elliptic PDE for h(6,(j)) on S using standard finite differencing methods. 



Rewriting (70) in the form 

G = G (/i, duh, duvh; -fij, Kij, dkjij) = , (73) 

the expansion G is evaluated on a trial surface, then iterated using a Newton- Raphson method 
to solve 3 ■ 6h = —0, where J is the Jacobian matrix. The drawback of this method is that 
it is not guaranteed to give the outermost marginally trapped surface. In practice however, 
this limitation can be easily overcome by either a single good initial guess, or multiple less 
accurate initial guesses. 

5.6.2. Masses and Momenta Two distinct measures of mass and momenta are available 
in relativistic spacetimes. First, ADM mass and angular momentum evaluated as either 
surface integrals at infinity or volume integrals over entire hypersurfaces give a measure of 
the total energy and angular momentum in the spacetime. The module ML_ADMQuantities 
of the McLachlan code |166j uses the latter method, creating gridfunctions containing the 
integrand of the volume integrals |167j : 



M -- 

IGtt 



— [ d^x e^'t' ( IGvrE + Ai.A'^ - -K"^] - e^R 



J^ =^^^J' jj""^ e'''l'(^Ah + '^xW,K-^x^Aendk'f''' + 8nx^Sk^ 



(74) 
(75) 



on which the user can use the reduction functions provided by Carpet to perform the volume 
integral. We note that ML_ADMQuantities inherits directly from the BSSN variables stored in 
McLachlan rather than strictly from the base modules. As the surface terms required when 
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converting a surface integral to a volume integral are neglected, this procedure assumes that 
the integrals of D^e'^ and e^'^eij''x^ A^k over the boundaries of the computational domain 
vanish. The ADM mass and angular momentum can also be calculated from the variables 
stored in the base modules using the Extract module, as surface integrals |136j 



on a specified spherical surface, preferably one far from the center of the domain since these 
quantities are only properly defined when calculated at infinity. 

There are also the quasi-local measures of mass and angular momentum, from any AHs 
found during the spacetime. Both AHFinderDirect and AHFinder output the corresponding 



The module QuasiLocalMeasures implements the calculation of mass and spin 
multipoles from the isolated and dynamical horizon formalism |168l 1169] , as well as a number 
of other proposed formula for quasilocal mass, linear momentum and angular momentum 
that have been advanced over the years |170] . Even though there are only a few rigorous 
proofs that establish the properties of these latter quantities, they have been demonstrated 
to be surprisingly helpful in numerical simulations (see, e.g., |171] ). and are therefore an 
indispensable tool in numerical relativity. QuasiLocalMeasures takes as input a horizon 
surface, or any other surface that the user specifies (like a large coordinate sphere) and can 
calculate useful quantities such as the Weyl or Ricci scalars or the three-volume element of 
the horizon world tube in addition to physical observables such as mass and momenta. 

Finally, the module HydroAnalysis additionally locates the (coordinate) center of mass 
as well as the point of maximum rest mass density of a matter field. 

5. 6.3. Gravitational Waves One of the main goals of numerical relativity to date is modeling 
gravitational waveforms that may be used in template generation to help analyze data from 
the various gravitational wave detectors around the globe. The Einstein Toolkit includes 
modules for extracting gravitational waves via either the Moncrief formalism of a perturbation 
on a Schwarzschild background or the calculation of the Weyl scalar \E'4. 

The module Extract uses the Moncrief formalism |172] to extract gauge-invariant 
wave functions Q^^ and given spherical surfaces of constant coordinate radius. The 
spatial metric is expressed as a perturbation on Schwarzschild and expanded into a tensor 
basis of the Regge- Wheeler harmonics |173j described by six standard Regge- Wheeler 



functions {cf"',cf"',ht^"',H+^'^,K+^"',G+^"'}. From these basis functions the gauge- 




(76) 



(77) 



mass derived from the area of the horizon tuh = \/ Aj (167r). 



invariant quantities: 




2(£ + 2)! 



c 



,x£m 
1 




s 



r 




The Einstein Toolkit 



31 



Q 



£m 




(79) 



are calculated, where S = 1 — 2M/r and A = {£ — + 2) + 6M/r. These functions then 
satisfy the wave equations: 

'£{£ + !) 6M 
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£m 



-S 

-S 



Q 



1 /72M3 12M(£- l)(£ + 2) 
^-l){l + 2) 



0) 
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3Af 



r2A 



where r* = r+2Mln(r/2M— 1). Since these functions describe the 4-metric as a perturbation 
on Schwarzschild, the spacetime must be approximately spherically symmetric for the output 
to be interpreted as first-order gauge-invariant waveforms. 

For more general spacetimes, the module WeylScal4 calculates the complex Weyl scalar 
^1^4 = Cap-ysnf^fh^n^fh^ , which is a projection of the Weyl tensor onto components of a null 
tetrad. WeylScal4 uses the fiducial tetrad |174] . written in 3+1 decomposed form as: 

r = i (u^ + (82) 



P) 



1 



^3) 



where is the unit normal to the hypersurface. The spatial vectors {f^, 6'^, 0^} are 
created by initializing f'^ = {0,x*}, 0^ = {0, —y,x,0}, and 9^ = {0, y/^Y^ek£m(p^r"^}, then 
orthonormalizing starting with 0* and invoking a Gram-Schmidt procedure at each step to 
ensure the continued orthonormality of this spatial triad. 

The Weyl scalar \E'4 is calculated explicitly in terms of projections of the 3-Riemann 
tensor onto a null tetrad, such that 



— f) i k — f\ 

m n-'n m j 



i5) 



For a suitably chosen tetrad, this scalar in the radiation zone is related to the strain of the 
gravitational waves since 



h = — ihs 



dt' 



^Adt" 



(86) 



While the waveforms generated by Extract are already decomposed on a convenient 
basis to separate modes, the complex quantity ^^4 is provided by WeylScal4 as a complex 
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grid function. For this quantity, and any other real or complex grid function, the module 
Multipole interpolates the field u{t,r,6,(f)) onto coordinate spheres of given radii and 
calculates the coefficients 

C'"" {t, r) = J ,Y;^u{t, r, 6, <j>ydn (87) 
of a projection onto spin-weighted spherical harmonics sY^m- 

5.6.4- Object tracking We provide a module (PunctureTracker) for tracking BH positions 
evolved with moving puncture techniques. It can be used with (CarpetTracker) to have the 
mesh refinement regions follow the BHs as they move across the grid. The BH position is 
stored as the centroid of a spherical surface (even though there is no surface) provided by 
SphericalSurf ace. 

Since the punctures only move due to the shift advection terms in the BSSN equations, 
the puncture location is evolved very simply as 

where is the puncture location and /3* is the shift. Since the puncture location usually 
does not coincide with grid points, the shift is interpolated to the location of the puncture. 



Equation ((88)) is implemented with a simple first-order Euler scheme, accurate enough for 
controlling the location of the mesh refinement hierarchy. 

Another class of objects which often needs to be tracked are neutron stars. Here is it 
usually sufficient to locate the position of the maximum density and adapt AMR resolution 
in these regions accordingly, coupled with the condition that this location can only move at 
a specifiable maximum speed. 

5.6.5. Other analysis modules The remaining analysis capabilities of the Einstein Toolkit 
span a variety of primarily vacuum-based functions. First, modules are provided to calculate 
the Hamiltonian and momentum constraints which are used to monitor how well the 
evolved spacetime satisfies the Einstein field equations. Two modules, ADMConstraints and 
ML_ADMConstraints provide these quantities. Both calculate these directly from variables 



stored in the base modules described in Sec. 5.1, explicitly written as: 

H = R- K'jKi, + -IQttE (89) 

Mi = VjKi^ - ViK - SnSi (90) 

where Si = (Tjo — P^Tij). The difference between these modules lies in how they access 
the stress energy tensor T^^., as the module ADMConstraints uses a deprecated functionality 
which does not require storage for T^y. 

Finally, ADMAnalysis calculates a variety of derived spacetime quantities that are often 
useful in post-processing such as the determinant of the 3-metric det7, the trace of the 
extrinsic curvature K, the 3-Ricci tensor in Cartesian coordinates TZij and its trace 7^, as 
well as the 3-metric and extrinsic curvature converted to spherical coordinates. 
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5. 7. Simulation Domain, Symmetries, Boundaries 

5.7.1. Domains and Coordinates. Cactus distinguishes between the physical domain, which 
hves in the continuum, and discrete domain, which consists of a discrete set of grid points. 
The physical domain is defined by its coordinate extent and is independent of the numerical 
resolution; in particular, the boundary of the physical domain has a width of zero (and is 
thus a set of measure zero). The discrete domain is defined indirectly via a discretization 
procedure that specifies the number of boundary points, their location with respect to the 
physical boundary, and either the grid spacing or the number of grid points spanning the 
domain. This defines the number and location of the grid points in the discrete domain. 
The discrete domain may have grid points outside of the physical domain, and may have a 
non-zero boundary width. This mechanism ensures that changes in the numerical resolution 
do not affect the extent of the physical domain, i.e., that the discrete domains converge to 
the physical domain in the limit of infinite resolution. The Einstein Toolkit provides the 
CoordBase thorn that facilitates the definition of the simulation domain independent of the 
actual evolution thorn used, allowing it to be specified at run time via parameters in the 
same way that parameters describing the physical system are specified. CoordBase exposes 
a public runtime interface that allows other thorns to query the domain description in a 
uniform way. This is used by Carpet to query CoordBase for the discrete grid when creating 
the hierarchy of grids, automatically ensuring a consistent grid description between the two 
thorns. Evolution thorns such as McLachlan use the domain information to decide which 
points are evolved and therefore require the evaluation of the right-hand-side expression, and 
which ones are set via boundary or symmetry conditions. 

5.7.2. Symmetries and Boundary Conditions. The Einstein Toolkit includes two thorns. 
Boundary and SymBase, to provide a generic interface to specify and implement boundary 
and symmetry conditions. The toolkit includes built-in support for a set of reflecting or 
rotating symmetry conditions that can be used to reduce the size of the simulation domain. 
These symmetries include periodicity in any of the coordinate directions (via the Periodic 
module), reflections across the coordinate planes (via the Reflection module), 90° and 180° 
rotational symmetries (via the RotatingSymmetryQO and RotatingSymmetrylSO modules 
respectively) about the z axis, and a continuous rotational symmetry (via the Cartoon2D 
thorn) |175j . Cartoon2D allows fully three dimensional codes to be used in axisymmetric 
problems by evolving a slice in the y = plane and using the rotational symmetry to 
populate ghost points off the plane (see Figure [s]) . 

In applying symmetries to populate ghost zones, the transformation properties of 
tensorial quantities (including tensor densities and non-tensors such as Christoffel symbols) 
are correctly taken into account, just as they are in the interpolation routines present in 
Cactus. Thus, symmetries are handled transparently from the point of view of user modules 
(see Figure |6] for an illustration). 

The Boundary thorn serves as a registry for available boundary conditions and provides 
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Figure 5. Grid layout of a simulation using Cartoon2D. The z-axis is the axis of rotational 
symmetry. Image courtesy of Denis PoUney. 




Figure 6. Iterative transformation of a point x in quadrant 3 to the corresponding point 
x" for which there is actual data stored. In this example, two reflection symmetries along 
the horizontal and vertical axis are present. Notice how the vector components change in 
transformations A and B. 
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basic scheduling to enforce all requested boundary conditions at the proper times. It also 
provides a basic set of boundary conditions to be used by user thorns. The "flat" boundary 
conditions often used for hydrodynamic variables that approach an atmosphere value fall in 
this category. More complicated boundary conditions are often implemented as modifications 
to the evolution equations and are not handled directly by Boundary. Examples are the 
radiative (Sommerfeld) and extrapolation boundary conditions provided by thorn NewRad. 

5.7.3. Adaptive Mesh Refinement The Einstein toolkit currently supports feature-based 
mesh refinement, which is based on extracting quantities such as the locations of BHs 
or NSs and then constructing a mesh hierarchy (stacks of refined regions) based on the 
locations, sizes, and speeds of these objects. This allows tracking objects as they move 
through the domain. One can also add or remove stacks if, for instance, the number of 
objects changes. Full AMR based on a local error estimate is supported by Carpet, but the 
Einstein Toolkit does not presently provide a suitable regridding thorn to create such a grid. 
If initial conditions are constructed outside of Carpet (which is often the case), then the 
initial mesh hierarchy has to be defined manually. In order to facilitate the description of the 
mesh hierarchy the Einstein toolkit provides two regridding modules in the CarpetRegrid 
and CarpetRegrid2 thorns. Both thorns primarily support box-in-box type refined meshes, 
which are well suited to current binary BH simulations in which the high-resolution regions 
are centered on the individual BHs. Figure [7] shows a typical set of nested boxes during the 
inspiral phase of a binary BH merger simulation. 




x=0 



Figure 7. Nested boxes following the individual BHs in binary BH merger simulation (see 
Section 6.2 1, with the location of the individual BHs found by PunctureTracker. The 



innermost three of the nine levels of mesh refinement used in this simulation are shown. 
Notice the use of RotatingSymmetrylSO to reduce the computational domain. 



CarpetRegrid provides a number of different ways to specify the refined regions, e.g., as 
a set of boxes centered around the origin or as an explicit list of regions that make up the grid 
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hierarchy. Traditionally, groups using CarpetRegrid have employed auxiliary thorns that are 
not part of the Einstein Toolkit to create this list of boxes based on information obtained 
from apparent horizon tracking or other means. CarpetRegrid2 provides a user-friendly 
interface to define sets of nested boxes that follow BHs or other tracked objects. Object 
coordinates are updated by CarpetTracker, which provides a simple interface to the object 
trackers PunctureTracker and NSTracker (see section 5.6.4) in order to have the refined 
region follow the moving objects. CarpetRegrid2 contains code to handle the vr-symmetry 
provided by RotatingSymmetrylSO, enforcing the symmetry on the resulting grid layout (see 
Figure [s]) . 
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Figure 8. Grid layout created by CarpetRegrid2. In tlris example we use one ghost point, 
one boundary point, and two buffer points as well as RotatingSymmetrylSO. There are 
two refinement levels present, a coarse one represented by big red circles and a fine one 
represented by small black circles. The filled black circle is the single point specified by 
the user. CarpetRegrid2 surrounded it with a layer of buffer points, indicated by the cyan 
filled circles. The open circles are ghost and boundary points which are maintained by 
Carpet. The presence of the 7r-symmetry forces CarpetRegrid2 to create the tiny region to 
the bottom left of the grid. It serves only as a source for the boundary condition. 



6. Examples 

To demonstrate the properties of the code and its capabilities, we have used it to simulate 
common astrophysical configurations of interest. Given the community-oriented direction of 
the project, the parameter files required to launch these simulations and a host of others 
are included and documented in the code releases, along with the data files produced by 
a representative set of simulation parameters to allow for code validation and confirmation 
of correct code performance on new platforms and architectures. As part of the internal 
validation process, nightly builds are checked against a set of benchmarks to ensure that 
consistent results are generated with the inclusion of all new commits to the code. 

The performance of the Toolkit for vacuum configurations is demonstrated through 
evolutions of single, rotating BHs and the merger of binary black hole configurations 



(sections |6.1 and 6.2, respectively). Linear oscillations about equilibrium for an isolated 
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NS are discussed in section 6.3 and the collapse of a NS to a BH, including dynamical 



formation of a horizon, in section 6.4 Finally, to show a less traditional application of the 



code, we show its ability to perform cosmological simulations by evolving a Kasner spacetime 



(see section 6.5 ) 



6.1. Spinning BH 

As a first example, we perform simulations of a single distorted rotating BH. We use 
TwoPunctures to set up initial data for a single puncture of mass Mbh = 1 and 
dimensionless spin parameter a = Si^^/M^^i = O-'''- Evolution of the data is performed 
by McLachlan, apparent horizon finding by AHFinderDirect and gravitational wave 
extraction by WeylScal4 and Multipole. Additional analysis of the horizons is done 
by QuasiLocalMeasures. The runs were performed with fixed mesh refinement provided 
by Carpet, using 8 levels of refinement on a quadrant grid (symmetries provided by 
Ref lectionSymmetry and RotatingSymmetrylSO). The outer boundaries were placed at 
R = 256M. We performed runs at 3 different resolutions: the low resolution was 
0.024M(3.072M), medium was 0.016M(2.048M) and high was 0.012M(1.536M), where the 
numbers refer to the resolution on the finest (coarsest) grid. The runs where performed using 
the tapering evolution scheme in Carpet to avoid interpolation in time during prolongation. 
The initial data corresponds to a rotating BH perturbed by a Brill wave and, as such, has a 
non-zero gravitational wave content. We evolved the BH using 4th-order finite differencing 
from T = OM until it had settled down to a stationary state at T = 120M. 

Figure [9] shows the £ = 2,m = mode of r\l/4 extracted at -R = 30M, and its numerical 
convergence. In the top plot the black (solid) curve is the real part and the blue (dashed) curve 
is the imaginary part of r\l/4 for the high resolution run. Curves from the lower resolution are 
indistinguishable from the high resolution curve at this scale. In the bottom plot the black 
(solid) curve shows the absolute value of the difference between the real part of the medium 
and low resolution waveforms while the blue (dashed) curve shows the absolute value of the 
difference between the high and medium resolution waveforms in a log-plot. The red (dotted) 
curve is the same as the blue (dashed) curve, except it is scaled for 4th order convergence. 
With the resolutions used here this factor is (0.016^ - 0.024^) / (0.012^ - 0.016^) ^ 5.94. 



Figure 10 shows similar plots for the £ = 4, m = mode of r\l/4, again extracted at 
R = 30M. The top plot in this case shows only the real part of the extracted waveform but 
for all three resolutions (black solid curve is high, blue dashed curve is medium and red dotted 
curve is low resolution). Since the amplitude of this mode is almost a factor of 20 smaller 
than the £ = 2,m = mode there are actually small differences visible between resolutions 
in the beginning of the waveform. The bottom plot shows the convergence of the real part 
of the £ = 4, m = mode (compare with the bottom plot in Figure [9]) and demonstrates 
that even though the amplitude is much smaller we still obtain close to perfect fourth-order 
convergence. 



In addition to the modes shown in Figure [9] and 10 we note that the extracted 




Figure 9. The extracted £ = 2,m — mode of ^4 as function of time from tiie higii 
resolution run (top plot). The extraction was done at i? = 30A/. Shown is both the real 
(solid black curve) and the imaginary (dashed blue curve) part of the waveform. At the 
bottom, we show the difference between the medium and low resolution runs (solid black 
curve), between the high and medium resolution runs (dashed blue curve), and the scaled 
difference (for 4th order convergence) between the medium and low resolution runs (dotted 
red curve) for the real part of the £ = 2, m = waveforms. 
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Figure 10. Real part of the extracted £ — 4:,ni — mode of ^'4 as function of time (top 
plot) for the high (solid black curve) , medium (dashed blue curve) and low (dotted red curve) 
resolution runs. The extraction was done at i? = 30Af. The bottom plot shows for the real 
part of the ^ = 4, m = waveforms the difference between the medium and low resolution 
runs (solid black curve) , the difference between the high and medium resolution runs (dashed 
blue curve) as well as the scaled (for 4th order convergence) difference between the medium 
and low resolution runs (dotted red curve). 
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£ = 4, m = 4 mode is non-zero due to truncation error, but shows fourth-order convergence 
to zero with resolution (this mode is not present in the initial data and is not excited during 
the evolution). Other modes are zero to round-off due to symmetries at all resolutions. 

Since there is non-trivial gravitational wave content in the initial data, the mass of the 



BH changes during its evolution. In figure 11 , we show in the top plot the irreducible mass as 
calculated by AHFinderDirect as a function of time at the high (black solid curve), medium 
(blue dashed curve) and low (red dotted curve) resolutions. The inset shows in more detail the 
differences between the different resolutions. The irreducible mass increases by about 0.3% 
during the first 40M of evolution and then remains constant (within numerical error) for the 
remainder of the evolution. The bottom plot shows the convergence of the irreducible mass 
by the difference between the medium and low resolutions (black solid curve), the difference 
between the high and medium resolutions (blue dashed curved) as well as the scaled difference 
between the high and medium resolutions for fourth-order (red dotted curve) and third-order 
(green dash-dotted curve). The convergence is almost perfectly fourth-order until T = 50M, 
then better than fourth-order until T = 60M, and finally between third-order and fourth- 
order for the remainder of the evolution. The lack of perfect fourth-order convergence at 
late times may be attributed to non-convergent errors from the puncture propagating to the 
horizon location at the lowest resolution. 



Finally, in Figure 12 we show the total mass (top plot) and the change in the spin, 
AS = S{t) — S{t = 0), as calculated by QuasiLocalMeasures. In both cases the black 
(solid) curve is for high, blue (dashed) for medium and red (dotted) for low resolution. Since 
the spacetime is axisymmetric the gravitational waves cannot radiate angular momentum. 
Thus any change in the spin must be due to numerical error and AS should converge to zero 



with increasing resolution. This is clearly shown in the bottom plot of Figure 12 ; the green 
(dash-dotted) curve (the high resolution result scaled by a factor of 1.78 for second-order 
convergence to the resolution of the medium resolution) and the blue (dashed) curve are on 
top of each other. Since the QuasiLocalMeasures thorn uses an algorithm which is only 
second-order accurate overall, this is the expected result. The increase of about 0.22% in the 
mass of the BH is caused solely by the increase in the irreducible mass. 

6.2. BH Binary 

To demonstrate the performance in the code for a current problem of wide scientific interest, 
we have evolved a non-spinning equal-mass BH binary system. The initial data represent a 
binary system in a quasi-circular orbit, with an initial separation chosen to be r = 6M so we 
may track the later inspiral, plunge, merger and ring down phases of the binary evolution. 
Table [T] provides more details about the initial binary parameters used to generate the initial 
data. The TwoPunctures module uses these initial parameters to solve ((To|), the elliptic 



Hamiltonian constraint for the regular component of the conformal factor (see section 5.2.3). 
The spectral solution for this example was determined by using [nA,nB,n^] = [28,28,14] 
collocation points, and, along with the Bowen-York analytic solution for the momentum 
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Figure 11. The top plot shows the irreducible mass of the apparent horizon as a function 
of time at low (black solid curve), medium (blue dashed curve) and high (red dotted curve) 
resolutions. The inset is a zoom in on the y-axis to more clearly show the differences between 
the resolutions. The bottom plot shows the convergence of the irreducible mass. The black 
(solid) curve shows the difference between the medium and low resolution results, the blue 
(dashed) curve shows the difference between the high and medium resolution results. The 
red (dotted) and green (dash-dotted) show the difference between the high and medium 
resolutions scaled according to fourth and third-order convergence respectively. 
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Figure 12. The top plot shows the total mass and the bottom plot shows the change in 
spin (i.e. AS — S{t) — S{t — 0) of the BH as a function of time. In both plots the black 
(solid) curve is for high, blue (dashed) for medium and red (dotted) for low resolution. In 
the bottom plot the green (dash-dotted) curve shows the high resolution result scaled for 
second-order convergence. The agreement with the medium resolution curve shows that the 
change in spin converges to zero as expected. 



constraints, represents constrained GR initial data {'jij, Kij}. The evolution is performed by 
the McLachlan module. 

The simulation domain spans the coordinate range 

[[a;min,a;max], bmin,Z/max], [2;min,2max]] = [[0, 120] , [- 120, 120] , [0, 120]] , whcrc wc havc taken 
advantage of both the equatorial symmetry (implemented by the Ref lectionSymmetry 
module) and the 180° rotational symmetry around the z-axis, which we apply at the x = 
plane using the RotatingSymmetrylSO module. Carpet provides a hierarchy of refined 
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0.13808 0.47656 0.984618 



Table 1. Initial data parameters for a non-spinning equal mass BH binary. The punctures 
are located on the a;-axis at positions xi and X2, with puncture bare mass parameters 
mi = m2 = m, and momenta ±p. 



grids centered at each puncture. Here, we used 7 levels of refinement, where the box edge 
coordinate lengths are given by [128,32,16,8,4,2] in units of the total binary mass, which 
is set to unity. Note that overlapping boxes are automatically redefined by Carpet into one 
unique region before the domain decomposition takes place. 

Figure 13 shows the two puncture tracks throughout all phases of the binary evolution, 
provided by the PunctureTracker module. In the same plot we have recorded the intersection 
of the apparent horizon 2-surface with the z = plane every time interval t = lOM during 
the evolution. A common horizon is first observed at t = 116M. These apparent horizons 
were found by the AHFinderDirect module and their radius and location information stored 
as a 2-surface with spherical topology by the SphericalSurf ace module. The irreducible 
mass and dimensionless spin of the merged BH were calculated by the QuasiLocalMeasures 
module, and were found to be 0.647M and -0.243M-2, respectively. 

Two modules are necessary to perform the waveform extraction. The first one, 
WeylScal4, calculates the Weyl scalar \l/4 in term of the metric components and its 
derivatives; these were computed to be 4-th order accurate in this example. The second 
module, Multipole, interpolates the Weyl scalars onto spheres with centers and radii 
specified by the user, and performs a spherical harmonic multipole mode decomposition. 
Figure 13 shows the real and imaginary parts of the (/ = 2, m = 2) mode for \l/4 extracted 
on a sphere centered at the origin at -Robs = 60M . The number of grid points on the sphere 
was set to be [nQ,n^] = [120,240], which yields an angular resolution of 2.6 x 10~^ radians, 
and an error of the same order, since the surface integrals were calculated by midpoint rule 
- a first order accurate method. 

In order to evaluate the convergence of the numerical solution, we ran five simulations 
with different resolutions, and focus our analysis on the convergence of the phase and 
amplitude of the Weyl scalar \l/4. The mesh spacings adopted for the coarser grid 
in the AMR hierarchy for these different runs were {/iiow, ^med, ^medh, ^Wgh, /'-higher} = 
{2.0M, 1.5M, 1.25M, l.OM, 0.75M}, respectively, while the finer grid spacings can be easily 
found by dividing them by 2'^ for the kth level of mesh refinement. For this example, we set 
{hL, hLd, ^Ldh, ^iigh, ^MgheJ = {3.125M, 2.344M, 1.953M, 1.563M, 1.172M} x 10"^ for the 
finest grid in the different AMR hierarchies, respectively. 

Here, we consider the phase 0(t) and the amplitude A{t) of the Weyl scalar \l/4 at 
-Robs = 60M. In order to take differences between the numerical values at two different 
grid resolutions, we use an 8-th order accurate Lagrange operator to interpolate the higher- 
accuracy finite difference solution into the immediately coarser grid. We have experimented 
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Figure 13. In the left panel, we plot the tracks corresponding to the evolution of two 
punctures initially located on the x-axis at x = ±3. The solid blue line represents puncture 
1, and the dashed red line puncture 2. The circular dotted green lines are the intersections of 
the apparent horizons with the z — plane plotted every lOAf during the binary evolution. 
A common horizon appears at t = 116M. In the right panel, we plot the real (solid blue 
line) and imaginary (dotted red line) parts of the {I — 2, to = 2) mode of the Weyl scalar ^'4 
as extracted at an observer radius of i?obs = 60M. 



with 4-th and 6-th order as well, to evaluate the level of noise these interpolations could 
potentially introduce, but did not observe any noticeable difference and we report here on 
results from the higher-order option. 

In Figure 14, we show the convergence of the amplitude and phase of the Weyl scalar by 
plotting the logarithm of the absolute value of the differences between two levels of resolution. 
The differences clearly converge to zero as the resolution is increased. We also indicate on 
both plots the time at which the gravitational wave frequency reaches u = 0.2/M. We 
follow a community standard, agreed to over the course of the NRAR pTB] collaboration, 
that constrains the numerical resolution so that the accumulated phase error is not larger 
than 0.05 radians at a gravitational wave frequency of u = 0.2/M. From the plot, we assert 
that the phase error between the higher and high resolutions and the one between high and 
medium-high resolutions satisfies this criterion, while the phase error between the medium- 
high and medium resolutions barely satisfies the criterion; and the one between medium and 
low resolutions does not. We conclude then that the three highest resolution runs do have 
sufficient resolution to extract waveforms for use in the construction of analytic waveform 
templates. 
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Figure 14. Weyl scalar amplitude (left panel) and phase (right panel) convergence. The 
long dashed red curves represent the difference between the medium and low-resolution runs. 
The short dashed orange curves show the difference between the medium-high and medium 
resolution runs. The dotted brown ones, the difference between high and medium-high 
resolutions, while the solid blue curves represent the difference between the higher and high 
resolution runs. The dotted vertical green line at t = 15AM indicates the point during the 
evolution at which the Weyl scalar frequency reaches u! = 0.2/ M. Observe that the three 
highest resolutions accumulate a phase error below the standard of 0.05 radians required by 
the NRAR collaboration. 



6.3. Linear oscillations of TOV stars 

The examples in the previous subsections did not include the evolution of matter within 
a relativistic spacetime. One interesting test of a coupled matter-spacetime evolution 
is to measure the eigenfrequencies of a stable TOV star (see, e.g., |177H181j ). These 
eigenfrequencies can be compared to values known from linear perturbation theory. 

We begin our simulations with a self-gravitating fluid sphere, described by a polytropic 
equation of state. This one-dimensional solution is obtained by the code described in 



section 5.2.5 , and is interpolated on the three-dimensional, computational evolution grid. 
This system is then evolved using the BSSN evolution system implemented in McLachlan 
and the hydrodynamics evolution system implemented in GRHydro. 

For the test described here, we set up a stable TOV star described by a polytropic 
equation of state p = Kp^ with K = 100 and F = 2, and an initial central density of 
Pc = 1.28 X 10^^. This model can be taken to represent a non-rotating NS with a mass of 
M = 1.4M0. The computational domain is a cube of length 640M with a base resolution 
of 2M (4M, 8M) in each dimension. Four additional grids refine the region around the star 
centered at the origin, each doubling the resolution, with sizes of 120M, 60M, 30M and 15M, 
resulting in a resolution of 0.125M (0.25M, 0.5M) across the entire star. 
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Figure 15. Evolution of the central density for the TOV star. Clearly visible is an initial 
spike, produced by the interpolation of the one-dimensional equilibrium solution onto the 
three-dimensional evolution grid. The remainder of the evolution however, the central density 
evolution is dominated by continuous excitations coming from the interaction of the stellar 
surface with the artificial atmosphere. 



In Figure 15 we show the evolution of the central density of the star over an evolution 
time of 1300M (6.5ms). The initial spike is due to the perturbation of the solution resulting 
from the interpolation onto the evolution grid. The remaining oscillations are mainly due 
to the interaction of the star and the artificial atmosphere and are present during the whole 
evolution. Given enough evolution time, the frequencies of these oscillations can be measured 
with satisfactory accuracy. 

In Figure 16 we show the power spectral density (PSD) of the central density 
oscillations computed from a full 3D relativistic hydrodynamics simulation, compared to 
the corresponding frequencies as obtained with perturbative techniques (kindly provided 
by Kentaro Takami and computed using the method described in |182j ). The PSD was 
computed using the entire time series of the high-resolution run, by removing the linear 
trend and averaging over Hanning windows overlapping half the signal length after padding 
the signal to five times its length. The agreement of the fundamental mode and first three 
overtone frequencies is clearly visible, but are limited beyond this by the finite numerical 
resolution. Higher overtones should be measurable with higher resolution, but at substantial 
computational cost. 

Within this test it is also interesting to study the convergence behavior of the coupled 
curvature and matter evolution code. One of the variables often used for this test is the 
Hamiltonian constraint violation. This violation vanishes for the continuum problem, but is 
non-zero and resolution-dependent in discrete simulations. The expected rate of convergence 
of the hydrodynamics code lies between 1 and 2. It cannot be higher than 2 due to the 



The Einstein Toolkit 



47 




0123456789 10 

/ [kHz] 



Figure 16. Eigenfrequency mode spectrum of a TOV star. Shown is the power spectral 
density of the central matter density, computed from a full 3D relativistic hydrodynamics 
simulation and compared to the values obtained by perturbation theory. The agreement of 
the frequencies of the fundamental mode and the first three overtones is clearly visible. 



directional flux-split algorithm which is of second order. Depending on solution itself, the 
hydrodynamics code is only of first order in particular regions, e.g., at extrema (like the 
center of the star), or at the stellar surface. 



Figure 17 shows the order of convergence of the Hamiltonian constraint violation, using 
the three highest-resolution runs, at the stellar center and a coordinate radius of r = 5M 
which is about half way between the center and the surface. The observed convergence rate 
for most of the simulation time lies between 1.4 and 1.5 at the center, and between 1.6 and 2 
at r = 5M, consistent with the expected data-dependent convergence order of the underlying 
hydrodynamics evolution scheme. 

6.4- Neutron star collapse 

The previous examples dealt either with preexisting BHs, either single or in a binary, or 
with a smooth singularity free spacetime, as in the case of the TOV star. The evolution 
codes in the toolkit are, however, also able to handle the dynamic formation of a singularity, 
that is follow a neutron star collapse into a BH. As a simple example of this process, we 



study the collapse of a non-rotating TOV star. We create initial data as in section 6.3 
using pc = 3.154 x 10~^ and Kid = 100, F = 2, yielding a star model of gravitational mass 
1.67 M0, that is at the onset of instability. As is common in such situations (e.g., [95j), we 
trigger collapse by reducing the pressure support after the initial data have been constructed 
by lowering the polytropic constant Kid from its initial value to K = 0.98 Kid = 98. To 
ensure that the pressure-depleted configuration remains a solution of the Einstein constraint 
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Figure 17. Convergence factor of Hamiltonian constraint violation at r = CM and r — 5M. 
The observed convergence order of about 1.5 at the center of the star is lower then the 
general second order of the hydrodynamics evolution scheme. This is expected because the 
scheme's convergence rate drops to first order at extrema or shocks, like the stellar center or 
surface. Consequently, the observed convergence order about half way between the stellar 
center and surface is higher than 1.5, but mostly below 2. 



equations (89) in the presence of matter, we rescale the rest mass density p such that the 
total energy density T„„ does not change: 

p' + K{p'y = p + K,,,p\ (91) 

Compared to the initial configuration, this rescaled star possesses a slightly higher central 
density and lower pressure. This change in K accelerates the onset of collapse that would 
otherwise rely on being triggered by numerical noise, which would not be guaranteed to 
converge to a unique solution with increasing resolution. In order to resolve the star as 
well as to push the outer boundary far enough away (so that the star and the numerical 
outer boundary are not in causal contact during the simulation) we employ a fixed mesh 
refinement scheme. The outermost box has a radius of Rq = 204.8 Mq and a resolution of 
3.2 Mq (2.4 Mq, 1.6 Mq, 0.6 Mq for higher convergence levels). Around the star, centered 
about the origin, we stack 5 extra boxes of approximate size 8x2^ Mq for < £ < 4, where 
the resolution on each level is twice that of the surrounding level. In order to resolve the 
large density gradients developing during the collapse, two more levels with radii 4 Mq and 
2Mq are placed inside the star. We use the PPM reconstruction method and the HLLE 
Riemann solver to obtain second-order convergent results in smooth regions. Due to the 
presence of the density maximum at the center of the star and the non-smooth atmosphere 
at the edge of the star, we expect the observed convergence rate to be somewhat lower than 
second order, but higher than first order. In Figure 18, we plot the approximate coordinate 
size of the star as well as the circumferential radius of the apparent horizon that eventually 
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Figure 18. Coordinate radius of the surface of the collapsing star and radius of the 
forming apparent horizon. The stellar surface is defined as the point where p is 100 
times the atmosphere density. R is the circumferential radius of the apparent horizon and 
Rg = 2 Mi, = 2 X 1.63 Mq. An apparent horizon forms at a time roughly equal to when 
the mass of the star is enclosed in its gravitational radius, forming a black hole and causally 
disconnecting the evolution in the interior from the outside spacetime. The lower x-axis 
displays time in code units where Mq = G = c = 1, and the upper x-axis shows the 
corresponding physical time using 1 Mq = 4.93 /zs. 




Figure 19. Convergence factor for the Hamiltonian constraint violation at the center of the 
collapsing star. We plot convergence factors computed using a set of 4 runs covering the 
diameter of the star with sa 60, 80, 120, and 240 grid points. The units of time on the upper 
and lower x-axes are identical to those of Figure [TSj 
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forms in the simulation. The apparent horizon is first found at approximately the time when 
the star's coordinate radius approaches its Schwarzschild radius, though one needs to keep 
in mind that the Schwarzschild radius is a circumferential radius, whereas the meaning of 
the coordinate radius in our BSSN calculation is likely somewhat different. In Figure 19, we 
display the convergence factor obtained from 



Hhi — — h2 



(92) 



Hh2 — /if — /I3 ' 

for the Hamiltonian constraint violation at the center of the collapsing star. Here Hh^ is the 



Hamiltonian constraint violation (89) at the center of the star for a run with resolution hi. 
Up to the time when the apparent horizon forms, the convergence order is an expected ~ 1.5. 
At later times, the singularity forming at the center of the collapsing star renders a pointwise 
measurement of the convergence factor at the center impossible. 

6. 5. Cosmology 

The Einstein Toolkit is not only designed to evolve compact-object spacetimes, but also to 
solve the initial-value problem for spacetimes with radically different topologies and global 
properties. In this section, we illustrate the evolution of an initial-data set representing a 
constant-t section of a spacetime from the Gowdy T^ class |183l I184j . Models in this class 
have the line element: 

ds^ = r-i/2eV2(_^^2 ^ ^^2) ^ rle^'idx + Qdyf + 6"^%'] (93) 

defined on a 3-torus —Xq < x < Xq, —yo < y < yo, —Zq < z < zq, with the functions P, 
Q and A to be determined by the Einstein equations. For P = Q = \ = 0, a coordinate 
transformation t = 4/3 r^/^ (plus a rescaling of the spatial coordinates) casts the line element 
into the form: 

ds^ = -dt^ + t^'\dx^ + dy^) + ^94) 

which represents the familiar Kasner spacetime for a homogeneous but anisotropically 
expanding universe. In the 34-1 decomposition described above, this reads: 

a{t) = 1 (95) 
P\t) =0 (96) 
7.,(t) =diag(t^/^t^/^^-^/=^) (97) 

K,,{t)= -diag(^^^/^^t^/^^r^/^) (98) 



In Figure 20, we show the full evolution of the t = 1 slice of spacetime (94), along with 
the associated error for a sequence of time resolutions. 





Figure 20. Top: the evolution of a vacuum spacetime of the type ( 93 ), with P 
the initial data are chosen as — 6ij and Kij = diag(— 2/3, 



= A = 0; 
^2/3,1/3). Bottom: the 
numerical error for a sequence of four time resolutions dt — [0.0125,0.025,0.05,0.1]; the 
errors are scaled according to the expectation for fourth-order convergence. 



7. Conclusion and Future Work 

In this article, we described the Einstein Toolkit, a collection of freely available and easy- 
to-use computational codes for numerical relativity and relativistic astrophysics. The code 
details and example results present in this article represent the state of the Einstein Toolkit 
in its release ET_2011_05 "Curie," released on April 21, 2011. 

The work presented here is but a snapshot of the Einstein Toolkit's ongoing development, 
whose ultimate goal it is to provide an open-source set of robust baseline codes to realistically 
and reproducibly model the whole spectrum of relativistic astrophysical phenomena 
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including, but not limited to, isolated black holes and neutron stars, binary black hole 
coalescence in vacuum and gaseous environments, double neutron star and neutron star 
- black hole mergers, core-collapse supernovae, and gamma-ray bursts. 

For this, much future work towards including proper treatments of magnetic fields, more 
complex equations of state, nuclear reactions, neutrinos, and photons will be necessary and 
will need to be matched by improvements in infrastructure (e.g., more flexible AMR on 
general grids) and computing hardware for the required fully coupled 3-D, multi-scale, multi- 
physics simulations to become reality. These tasks, as well as the others mentioned below, 
are likely to occupy a great deal of the effort spent developing future versions of the Einstein 
Toolkit over the next few years. 

Without a doubt, collapsing stars and merging BH-NS and NS-NS binaries must be 
simulated with GRMHD to capture the effects of magnetic fields that in many cases will alter 
the simulation outcome on a qualitative level and may be the driving mechanisms behind 
much of the observable EM signature from GRBs (e.g., |185j ) and magneto-rotationally 
exploding core-collapse supernovae (e.g., |186] ). To date, all simulations that have taken 
magnetic fields into account are still limited to the ideal MHD approximation, which 
assumes perfect conductivity. Non-ideal GRMHD schemes are just becoming available (see, 
e.g., |187l 1188] ). but have yet to be implemented widely in many branches of numerical 
relativity. 

Most presently published 3D GR(M)HD simulations, with the exception of recent work 
on massive star collapse (see, e.g., [87]) and binary mergers (see, e.g., |18]), relied on simple 
zero-temperature descriptions of NS stellar structure, with many assuming simple polytropic 
forms. Such EOSs are computationally efficient, but are not necessarily a good description 
for matter in relativistic astrophysical systems. The inclusion of finite-temperature EOSs, 
derived from the microphysical descriptions of high-density matter, will lead to qualitatively 
different and much more astrophysically reliable results (see, e.g., [87]). In addition, most 
GR(M)HD studies neglect transport of neutrinos and photons and their interactions with 
matter. Neutrinos in particular play a crucial role in core-collapse supernovae and in 
the cooling of NS-NS merger remnants, thus they must not be left out when attempting 
to accurately model such events. Few studies have incorporated neutrino and/or photon 
transport and interactions in approximate ways (see, e.g., [HI |66l ETl 1189] ). 

Besides new additions of physics modules, existing techniques require improvement. One 
example is the need for the gauge invariant extraction of gravitational waves from simulation 
spacetimes as realized by the Cauchy Characteristic Extraction (CCE) technique recently 
studied in |112[ I190[ 1191] . The authors of one such CCE code |190j have agreed to make 
their work available to the whole community by integrating their CCE routines into the 
Einstein Toolkit release 2011_11 "Maxwell," which will be described elsewhere. 

A second much needed improvement of our existing methods is a transition to 
cell-centered AMR for GR hydrodynamic simulations, which would allow for exact flux 
conservation across AMR interfaces via a refiuxing step that adjusts coarse and/or fine grid 
fluxes for consistency (e.g., (ITT] ). This is also a prerequisite for the constrained transport 
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method |192] for ensuring the divergence-free condition for the magnetic field in a future 
implementation of GRMHD within the Einstein Toolkit. Work towards cell-centered AMR, 
refluxing, and GRMHD is underway and will be reported in a future publication. 

While AMR can increase resolution near regions of interest within the computational 
domain, it does not increase the convergence order of the underlying numerical methods. 
Simulations of BHs can easily make use of high-order numerical methods, with eighth-order 
convergence common at present. However, most GRMHD schemes, though they implement 
high-resolution shock-capturing methods, are limited to 2nd-order numerical accuracy in the 
hydrodynamic/MHD sector while performing curvature evolution with 4th-order accuracy or 
more. Higher order GRMHD schemes are used in fixed-background simulations (e.g., |193] ). 
but still await implementation in fully dynamical simulations. 

Yet another important goal is to increase the scalability of the Carpet AMR 
infrastructure. As we have shown, good scaling is limited to only a few thousand processes 
for some of the most widely used simulation scenarios. Work is in progress to eliminate 
this bottleneck |194| . On the other hand, a production simulation is typically composed of a 
large number of components, and even analysis and I/O routines have to scale well to achieve 
overall good performance. This is a highly non-trivial problem, since most Einstein Toolkit 
physics module authors are neither computer scientists nor have they had extensive training in 
parallel development and profiling techniques. Close collaboration with experts in these topics 
has been fruitful in the past and will be absolutely necessary for the optimization of Einstein 
Toolkit codes for execution on the upcoming generation of true petascale supercomputers on 
which typical compute jobs are expected to be running on 100,000 and more compute cores. 
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